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Abstract 

The  constitutive  relations  found  in  traditional  Navier-Stokes-based  computa¬ 
tional  fluid  dynamics  solvers  are  known  to  be  limited  in  altitude.  The  presence  of 
nonequilibrium  phenomena  beyond  what  these  methods  are  able  to  predict  becomes 
more  prevalent  at  higher  altitudes,  or  increasing  Knudsen  number.  The  bulk  viscosity, 
normally  assumed  to  be  zero  in  most  computational  fluid  dynamics  applications,  is 
examined  as  a  means  of  increasing  the  range  of  applicability  of  computational  fluid 
dynamics.  The  bulk  viscosity  model  used  was  from  recent  calculations  available  in 
the  literature,  from  a  new  anisotropic  potential  energy  surface,  and  is  restricted  to 
temperatures  below  2000  K.  The  normal  shock  problem  was  solved  for  Mach  numbers 
up  to  ten,  using  the  bulk  viscosity  model.  The  bulk  viscosity  provided  improvement 
in  the  agreement  with  observations  of  normal  shock  thickness  for  Mach  numbers  up  to 
ten.  Two  axisymmetric,  experimentally  observed  flows  were  solved  with  and  without 
the  bulk  viscosity,  and  compared  to  DSMC  solutions  of  a  previous  work.  Improvement 
of  surface  heat  transfer  agreement  with  observation  was  found  for  a  hollow-cylinder 
flare  axisymmetric  body.  Improvement  of  separation  point  prediction  was  found  for 
an  axisymmetric  double  cone. 
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Extending  CFD  Modeling  to  Near-Continuum  Flows 
Using  Enhanged  Thermophysigal  Modeling 

I.  Introduction 

The  concepts  of  maneuver  and  mass  have  been  redefined  by  the  sheer  speed  air 
and  space  vehicles  have  brought  to  warfare.  In  this  area,  major  technological 
battles  have  played  out  for  speed  and  range.  The  advantage  afforded  by  speed  in 
the  battlespace  is  well  known  [3:9-11;  4:8].  Air  Force  Chief  of  Staff,  General  Michael 
Moseley,  puts  speed  in  the  context  of  the  greater  mission  of  the  service: 

The  Air  Force  exists  to  fly,  hght  and  win — to  achieve  strategic,  opera¬ 
tional  and  tactical  objectives — unhindered  by  time,  distance  or  geography. 

. . .  We  provide  Global  Vigilance,  Global  Reach,  and  Global  Power. . . 

•  Global  Vigilance  is  the  persistent,  world-wide  capability  to  keep 
an  unblinking  eye  on  any  entity. . . 

•  Global  Reach  is  the  ability  to  move,  supply,  or  position  assets — 
with  unrivaled  velocity  and  precision — anywhere  on  the  planet 
•  Global  Power  is  the  ability  to  hold  at  risk  or  strike  any  target, 
anywhere  in  the  world,  and  project  swift,  decisive,  precise  effects. 

15: 1| 

Hypersonic  vehicles,  which  move  through  the  atmosphere  much  faster  than  the 
speed  of  sound,  have  provided  Global  Vigilance  through  space  launch,  and  Global 
Power  through  intercontinental  ballistic  missiles.  Both  these  applications  have  very 
high-speed,  short- duration  flights,  and  depend  on  rocket  propulsion.  As  successful  as 
these  types  of  vehicle  are,  there  has  long  been  a  quest  for  air-breathing  hypersonic 
vehicles.  As  imagined  by  designers,  these  vehicles  could  be  operated  from  runways  or 
underwing  stores  stations,  not  launchpads.  Global  power,  in  providing  more  precise 
effects  than  the  strategic  counterpart  mentioned  above,  may  become  the  domain  of 
these  new  vehicles. 

Global  vigilance  may  be  augmented  by  hypersonic  vehicles.  With  peer-competitors 
flexing  their  muscles  in  space,  the  Ghief  of  Staff  has  this  to  say:  “We  need  to  deploy 
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high- altitude,  high-speed,  air-breathing  systems  to  mitigate  risks  to  space-based  ca¬ 
pabilities.”  [5:8].  The  high- altitude,  high-speed  flight  regime  is  the  focus  of  this  work. 

In  producing  hypersonic  vehicles,  design  activities,  computer  simulations  and 
analyses  are  performed.  Ground  tests  are  conducted,  matching  as  many  of  the  pro¬ 
posed  flight  conditions  as  possible.  Actual  flight  research,  when  budgeted,  is  per¬ 
formed  to  validate  analytical  and  computer  models  and  reduce  their  uncertainty. 

The  importance  of  research  to  develop  this  technology,  and  more  importantly 
to  understand  the  physical  phenomena  associated  with  flying  at  high  speeds  and  high 
altitudes  cannot  be  stressed  enough.  The  United  States  Air  Force  (USAF)  must  be 
able  to  honestly  assess  the  effectiveness,  cost  and  risk  of  proposed  investments  in 
hypersonic  vehicles  and  weapons.  Through  these  research  programs,  data  for  testing 
flow  models  are  obtained,  and  better  models  can  be  bnilt.  Better  models  provide 
more  accnrate  predictions.  More  accnrate  predictions  can  reduce  risk,  uncertainty, 
and  cost.  Models  which  agree  with  a  particular  set  of  experiments  or  analysis  should 
not  be  treated  as  finished  works.  They  should  instead  be  snbject  to  refinement,  and 
re-validation,  especially  in  cases  where  vehicles  begin  to  be  used  beyond  their  original 
operating  envelope,  or  beyond  their  expected  life. 

Now  is  an  exciting  time  for  hypersonics  research  in  the  USAF.  The  X-51,  a  ve¬ 
hicle  designed  to  demonstrate  the  feasibility  of  snpersonic  combustion,  a  hypersonics- 
enabling  technology,  is  being  prepared  for  flight  test.  The  USAF  is  leading  the  world 
in  hypersonics  research.  In  fact,  at  the  46th  Aerospace  Sciences  Meeting  of  the  Amer¬ 
ican  Institnte  of  Aeronantics  and  Astronantics  (AIAA),  eight  hypersonics-related  ses¬ 
sions  were  chaired  by  Air  Force  personnel  [6].  A  joint  Australian-USAF  program,  the 
Hypersonics  International  Flight  Research  Experimentation  (HiFIRE),  is  currently  in 
gronnd  testing,  and  will  be  flown  within  the  year.  Gronnd  testing  done  at  National  Air 
and  Space  Administration  (NASA)  Langley  Research  Genter  (LaRG)  and  the  Galspan 
University  of  Buffalo  Research  Genter  (GUBRG),  combined  with  computational  flnid 
dynamics,  is  being  used  to  determine  the  geometry  of  the  flight  test  vehicle  [7]. 
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1.1  Continuum  Assumption  and  Altitude 

Most  aerodynamicists  are  used  to  thinking  of  air  as  a  continuous  medium,  in¬ 
finitely  divisible  into  smaller  pieces  with  identical  behavior.  The  reality  is,  in  fact, 
much  more  complicated.  Air  is  made  up  of  N2,  O2,  and  smaller  amounts  of  CO2,  Ar, 
and  other  ‘trace’  components.  Transfer  of  mass,  momentum,  and  energy  in  various 
forms  happen  on  a  molecular,  and  therefore  discrete,  level.  Continuum  representa¬ 
tions  are  built  from  the  microscopic,  making  assumptions  and  simplihcations  about 
the  collisional  processes  and  molecular  energy  transfer. 

A  “rough  guide”  [2:665]  for  the  dehnitions  of  the  flight  regimes  discussed  above 
can  be  found  in  Table  1 

Table  1:  Rarehed  Flow  Regimes  [2:665] 

Flight  regime  Range 

collisionless,  near-collisionless:  10.0  <  Kn 
transition:  0.1  <iFn<10.0 

near-hydrodynamic:  O.Ol^A'n^  0.1 

hydrodynamic:  Kn<  0.01 

By  this  guide,  using  body  dimension  as  the  basis  for  Kn,  the  present  work 
will  be  conhned  to  the  hydrodynamic  regime.  Thus,  this  work,  being  a  small  step 
toward  the  near-continuum  regime  will  start  from  within  the  continuum  regime.  The 
dehnition  of  high  Knudsen  number  will  be  Kn  >  0.001.  The  reason  for  this  dehnition 
will  be  explained  in  more  detail  in  Chapter  II.  Figure  1  depicts  mean  free  path,  as 
calculated  from  the  1976  US  Standard  Atmosphere  as  a  function  of  altitude.  With 
this  hgure,  and  a  local  length  scale,  a  Knudsen  number  can  be  determined. 

1.2  Near- Continuum  Flow  Considerations 

The  continuum  model,  traditionally  used  from  the  earth’s  surface  to  around 
80km,  is  based  on  the  familiar  continuum  formulations  of  conservation  of  mass,  mo¬ 
mentum  and  energy,  hereinafter  referred  to  collectively  as  the  Navier-Stokes  equations. 


3 


Figure  1:  Mean  free  path  from  1976  US  Standard  Atmosphere  [8] 
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The  higher  altitude,  more  rarehed  regime  has  been  the  domain  of  the  Boltzmann 
equation.  One  method  very  popular  today  of  solving  this  equation  is  the  Direct 
Simulation  Monte  Carlo  (DSMC)  method  due  to  Bird  [9].  The  main  problem  keeping 
this  method  from  being  applied  to  lower  altitude  problems  is  one  of  computational 
expense. 

Recently,  there  has  been  progress  in  bnilding  hybrid  methods  [10],  which  solve 
the  Navier-Stokes  eqnations  throngh  the  flowheld,  and  by  a  “continnum  breakdown” 
parameter  estimating  local  Knndsen  nnmber. 


Kri=^  (1.1) 

a  region  of  noneqnilibrinm  flow  is  identihed,  and  that  region  is  solved  by  a  DSMC 
solver,  and  the  two  are  conpled.  This  characteristic  length, L,  is  often  based  on  density 
gradient: 


Knp 


A'  dp 

p  dxi 


(1.2) 


Researchers  at  the  Air  Force  Research  Laboratory  have  shown  that  entropy 
generation  can  be  used  as  a  continuum  breakdown  parameter  [11]  for  assessing  the 
degree  of  non-eqnilibrium  present  in  a  flow. 

Several  researchers  [12]  [13]  [14]  have  taken  a  different  approach,  developing  gen¬ 
eralized  hydrodynamic  theories  which  go  beyond  the  Navier-Stokes  eqnations.  Most 
hydrodynamic  models  are  developed  from  the  Boltzmann  equation  by  using  some 
assnmptions,  simplihcations,  and  series  expansions  on  the  Maxwell-Boltzmann  distri- 
bntion  of  particles  in  velocity-space.  Chapman  and  Enskog  [15]  separately  developed 
similar  expansions  showing  that  a  zero-order  expansion  generates  the  Euler  equations, 
while  a  first-order  expansion  generates  the  Navier-Stokes  equations. 

The  natural  hrst  step  in  extending  these  equations  has  been  to  inclnde  more 
terms  in  the  expansion,  (e.g.  Burnett,  Super-Bnrnett,  BGK-Burnett  equations).  The 
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problem  with  these  methods  is  that  most  of  them  have  been  shown  to  violate  the 
second  law  of  thermodynamics  in  some  areas  [16].  The  goal  of  this  work  is  to  pro¬ 
vide  designers  with  the  simplest  possible  method  or  model  for  extending  the  current 
continuum  computational  fluid  dynamics  (CFD)  solvers  into  the  transition  regime. 

An  additional  consideration  for  high  Knudsen  number  flows  is  the  breakdown 
of  the  no-slip  condition.  This  condition  specifies  u  =  v  =  0  at  the  wall.  As  the 
gas  becomes  more  rarehed,  the  velocity  prohle  can  be  discontinuous,  and  the  gas 
temperature  can  be  different  from  the  adjacent  wall  temperature. 

1.3  Why  DSMC? 

As  mentioned  in  [17],  [18],  flows  with  appropriately  defined  local  Knudsen  num¬ 
bers  greater  than  0.3  cannot  be  appropriately  modeled  with  the  traditional  Navier- 
Stokes  approach.  The  continuum  constitutive  relations  and  their  usual  assumption  of 
isotropic  stresses  do  not  adequately  describe  the  physics  of  the  flow. 

Basing  Kn  on  the  overall  length  of  the  body  is  not  really  the  right  approach, 
as  the  mean  free  path  varies  throughout  hypersonic  flowhelds.  Bird  [9]  dehnes  L  for 
local  Knudsen  number  as 


dp/ dx 

The  traditional  continuum  description  of  flow  is  based  mostly  on  the  idea  that  at 
standard  conditions,  there  are  enough  molecules  present  that  the  flow  can  be  treated 
in  terms  of  properties  averaged  over  the  molecules,  as  conforming  to  smooth,  well- 
behaved  distribution  functions.  The  “macroscopic”  properties  at  a  “point,”  such 
as  density,  pressure,  temperature,  velocity,  et  cetera  are  actually  spatially  and/or 
time-averaged  properties,  over  a  “differential”  volume.  According  to  Bird,  when  con¬ 
structing  gradients  in  this  description,  what  is  actually  required  is  a  hnite  change  in 
that  property  over  several  mean  free  paths.  When  this  is  not  the  case,  the  standard 
continuum  model  breaks  down  [9:2]. 
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Bird  [9:4],  suggests  a  continuum  breakdown  parameter  (steady  flow  version): 


P 


A'  dp 

- s —  — 

2  p  dx 


(1.4) 


where  A'  is  mean  free  path,  s  is  a  ‘speed  ratio’  of  velocity  to  the  most  probable  speed, 
(2i?T)i/2.  Bird  sets  P  =  0.2  as  a  value  at  which  initial  breakdown  occurs. 

The  problem  of  predicting  breakdown  in  near-continuum  flow  is  still  of  interest 
today,  at  AFRL  [11].  Even  DSMC  has  low  Knudsen-number  limit:  a  point  at  which 
there  are  significant  fluctuations  in  the  statistics  of  the  algorithm  [9:21]. 


1.4  Molecular  Description  of  Flows 

The  molecular  description  of  flow  has  long  been  studied  as  the  iV-body  problem. 
The  major  difficulty  is  the  fact  that  iV  is  a  large  number,  and  modeling  forces  of  all 
the  bodies  on  each  other  is  extremely  computationally  expensive. 

Hamiltonian  dynamics  gives  very  interesting  results  for  these  types  of  prob¬ 
lems,  however  typically  very  few  molecules  are  used  in  any  analytical  work.  Many 
have  attempted  analytical  solutions  to  the  problem  for  gases  in  a  flow  of  interest  to 
aerodynamicists,  to  no  avail. 

Simulations  have  provided  some  insight,  especially  in  the  area  of  molecular 
dynamics,  see  [19],  for  example.  The  problem  again  being,  that  except  for  flows  with 
very  few  molecules  or  those  treated  without  wall  interactions,  the  vast  number  of 
bodies  needed  for  simulation  makes  problems  of  interest  in  hypersonics  intractable. 


1.5  Boltzmann’s  Equation 

Boltzmann’s  ideas  prompted  the  derivation  of  an  integro-differential  equation, 
now  referred  to  as  Boltzmann’s  Equation,  here  for  a  dilute  gas: 
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(9  (9  d 

—  {nf)  +  c-—{nf)  +  ¥-—{nf)  =  j  J  fl  -  f  fi)CradVLdci  (1.5) 

where  starred  quantities  are  post  collision,  n  is  the  number  density,  c  is  the  velocity 
class  of  one  collision  partner,  Ci  is  the  velocity  class  of  the  other  (this  is  a  two-body 
collision),  /  is  the  distribution  of  molecules  over  velocity  space,  associated  with  c,  /i 
is  the  same,  associated  with  Ci,  a  is  the  cross-section,  and  dVL  is  a  differential  solid 
angle.  F  is  the  net  effect  of  external  forces,  such  as  an  externally  applied  electric  held. 

The  term  on  the  right  hand  side  of  (1.5)  is  called  the  collision  integral,  itself  a 
very  difficult  quantity  to  compute.  Some  researchers  [20]  have  computed  this  integral 
over  actual  velocity  space  using  Monte  Carlo  integration,  and  treated  the  left  hand 
side  somewhat  similarly  to  modern  CFD  approaches,  to  study  normal  shock  structure. 

The  Boltzmann  equation  relates  collisions  sampled  from  continuous  velocity 
distribution  functions  to  time  and  spatial  rates  of  change  of  those  continuous  velocity 
distribution  functions.  One  could  imagine  then,  that  there  may  be  a  low-density 
limit  on  the  usefulness  of  the  Boltzmann  equation,  collisionless  Boltzmann  equation, 
or  DSMC,  in  terms  of  a  situation  where  there  are  not  enough  molecules  present  in 
the  volume  of  interest  to  compute  a  continuous  velocity  distribution  which  had  any 
meaning.  There  is  also  an  idea  of  ergodicity  inherent  in  this  equation.  That  is  to  say, 
spatial  and  time  derivatives  of  averaged  quantities  are  equated. 

1.6  Processes  of  DSMC 

The  basic  idea  in  DSMC  is  not  to  solve  the  Boltzmann  Equation  directly  by 
an  analytic  process,  but  rather  to  use  a  stochastic  model  to  simulate  the  Boltzmann 
Equation,  with  statistics  and  collision  processes  which  try  to  match  kinetic  theory. 

The  hrst  operation  in  DSMC  is  the  discretization  of  the  geometry  of  the  flowheld 
into  cells  and  assignment  of  boundary  conditions  just  as  in  CFD.  The  next  step  is 
gathering  information  about  the  flow  constituents,  such  as  molecular  weight,  rotation 


and  vibration  characteristic  temperatnres,  moments  of  inertia,  dissociation  energies, 
and  other  parameters  pertaining  to  the  flow  of  interest.  The  third  step  is  initializing 
the  flow,  usnally  to  some  equilibrinm  distribntion  abont  the  mean  freestream  velocity, 
using  a  pseudorandom  number  generator.  Pseudorandom  number  generators  will  be 
discussed  further  in  the  next  section. 

Rather  than  integrating  trajectories  deterministically,  the  main  algorithm  ap¬ 
plies  a  stochastic  model  to  determine  how  many  collisions  occur  within  each  cell, 
how  much  energy  is  transferred  to  the  various  modes  (translational,  rotational,  vibra¬ 
tional),  and  which  molecules  will  move  to  adjacent  cells.  Of  course,  this  model  should 
be  subject  to  constraints  of  conservation  of  particles,  momenta,  and  energy.  Critics 
of  DSMC  point  out  that  it  is  very  difficult  to  prove  within  the  simulation  that  these 
quantities  are  actually  conserved. 

Molecules  (or  atoms)  in  wall-bounded  cells  are  treated  according  to  an  acco¬ 
modation  coefficient,  with  one  extreme  being  diffuse  reflection  (somewhat  akin  to  a 
no-slip  condition)  and  the  other  extreme  being  specular  reflection.  Most  simulations 
discussed  in  [9]  are  accomplished  with  pure  diffuse  reflection.  Diffuse  reflection  refers 
to  the  situation  of  the  outgoing  velocity  of  the  collision  being  assigned  a  Maxwellian 
(or  Rayleigh  in  the  normal  direction)  distribution  referenced  to  the  wall,  whereas 
specular  reflection  refers  to  the  molecule  being  reflected  in  the  normal  direction  with 
the  same  speed,  and  retaining  its  tangential  speed.  This  model,  though  it  may  be 
beyond  the  no-slip  treatment  afforded  by  many  CFD  solvers,  is  quite  a  simplihcation 
of  what  may  actually  be  happening  at  the  wall,  in  terms  of  force  potential  and  energy 
transfer.  However,  Bird  [9]  gives  several  examples  where  this  is  reasonable. 

During  this  process,  the  program  systematically  takes  samples  of  particle  veloc¬ 
ity  and  other  parameters  at  user-dehned  intervals.  The  idea  is,  once  the  transients 
due  to  the  initial  conditions  have  died,  the  velocities  (and  time-averaged  quantities 
based  on  them)  reach  some  steady-state  distribution. 
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From  the  law  of  large  numbers,  if  some  quantity  obeys  a  particular  distribution 
over  the  entire  population,  the  cumulative  distribution  of  samples  approaches  that  of 
the  population  as  the  inverse  root  of  the  number  of  samples: 

liniM^largel  <X>  -E[<  X  >]\  ^  ^  (1.6) 

where  M  is  the  number  of  observations,  <  a;  >  is  the  average  over  the  samples,  and 
E[<  X  >]  is  the  expected  value  over  many  samples. 

One  can  imagine  this  process  getting  computationally  expensive.  One  thing 
DSMC  does  to  reduce  this  problem,  is  to  simulate  “pseudoparticles,”  which  represent 
many  actual  particles  (molecules,  atoms,  etc.).  There  are  assumptions  which  go  along 
with  this  idea,  the  most  common  being  that  all  the  particles  within  each  pseudoparticle 
are  in  equilibrium  with  each  other.  Bird  gives  data  to  support  the  idea  that  this  works, 
in  many  cases. 

1 . 7  Numerics 

DSMC  is  highly  dependent  on  a  “good”  pseudorandom  number  generator.  Gen¬ 
erators  which  repeat  too  often,  and/or  have  “clumping”  of  the  distributions  of  pseu¬ 
dorandom  numbers  can  effect  the  results,  or  greatly  increase  the  time  required  to 
obtain  a  solution  with  stationarity  of  the  averages  over  the  samples.  Some  workers 
refer  to  this  behavior  as  convergence,  but  it  bears  no  mathematical  resemblance  to 
convergence  of  an  algorithm  for  solving  a  deterministic  equation.  Optimally,  the  grid 
should  be  such  that  the  cell  size  is  about  one  third  of  the  mean  free  path. 

In  contrast  to  DSMC,  in  CFD,  the  solution  is  based  on  an  iterative  process.  Cell 
size  is  related  to  the  magnitude  of  gradients  to  be  captured.  Another  big  issue  with 
CFD  is  measuring  the  numerical  dissipation  renders  the  converged  solution  different 
from  the  actual  flow  being  modeled,  where  this  dissipation  is  required  by  many  CFD 
schemes  for  numerical  stability.  One  can  also  speak  of  grid  convergence,  meaning 
the  grid  resolution  at  which  the  solution  does  not  improve  with  grid  rehnement.  In 
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the  hypersonic  situation,  cell  sizes  have  to  make  sense  in  terms  of  an  appropriate 
breakdown  parameter,  as  mentioned  above.  If  the  cells  are  smaller  than  a  few  mean 
free  paths,  the  solution  can  become  less  physically  meaningful  [9:4], 

It  is  the  argument  of  Bird,  developer  of  DSMC,  that  it  is  not  physically  mean¬ 
ingful  to  compute  gradients  in  the  continuum  description  in  this  case  [9:4],  There  are 
other  arguments,  e.g.,  Lumpkin  [21:4],  Carr  [11:22-23]  which  submit  that  the  consti¬ 
tutive  model  could  be  redehned  for  the  continuum  transition  regime,  and  continuum 
solutions  could  be  obtained  from  these.  There  is  a  held  of  study  devoted  to  this, 
called  generalized  hydrodynamics  [14],  [12]. 

1.8  Computational  Cost 

Computational  cost  is  time  spent  per  solution,  either  in  total  CPU  time  ,  or 
“wall  time.”  Of  the  three  major  computational  methods  mentioned,  CFD  is  generally 
fastest.  The  molecular  dynamics  method  is  generally  considered  the  slowest.  One 
reason  for  this  is  the  total  number  of  computations  performed  per  time  step,  which 
can  be  related  to  the  number  or  grid  points,  and  in  the  particle  simulations,  to  the 
number  of  particles.  The  other  consideration  is  the  choice  of  a  time  step.  Since 
molecular  dynamics  is  a  deterministic,  mechanics-based  solver,  the  time  step  must  be 
small  enough  to  allow  integration  of  molecular  trajectories.  In  DSMC,  the  time  step 
is  not  as  restrictive,  but  still  must  be  less  than  the  mean  time  between  collisions.  The 
movement  is  broken  into  two  main  steps,  one  for  movement  of  particles,  and  a  second 
for  collision  events,  based  on  the  distribution  of  particles  at  the  end  of  the  movement 
step.  In  explicit-integration  CFD,  the  time  step  is  based  on  the  speed  a  wave  can 
travel  in  the  fluid  and  cell  size.  Information  propagates  through  the  grid  through  lo¬ 
cal  calculation  of  the  convective  derivative.  In  implicit-integration  CFD,  the  timestep 
can  be  theoretically  be  many  times  that  of  the  explicit.  Implicit  integration  solves 
the  time  derivative  globally,  allowing  information  to  propagate  much  farther  across 
the  grid.  The  global  time  derivative  can  be  thought  of  as  an  operator  or  matrix  to  be 
inverted  to  solve  for  the  values  in  the  grid  at  the  next  time  step.  The  DSMC  method 
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is  more  computationally  intensive  than  CFD,  however,  it  can  potentially  provide  a 
more  accurate  answer  in  regions  where  the  hydrodynamic  equations  being  solved  are 
not  valid.  Here  it  is  assumed  the  DSMC  software  is  given  the  input  paramters  and 
boundary  conditions  to  most  closely  match  the  physical  observations.  As  computers 
become  more  powerful,  the  higher  computational  expense  of  DSMC  over  CFD  may 
make  less  difference  to  designers.  Bird  writes  in  his  1994  book  [9]  that  all  his  simu¬ 
lations  took  less  than  24  hours.  Realistic  flows  of  interest  in  hypersonic  engineering 
are  still  more  than  one  day  per  solution  today. 

Alternative  solutions  from  the  CFD  world  are  beginning  to  show  some  success, 
however  [22],  [16].  There  may  be  more  overlap  in  the  “transition  region,”  meaning 
transition  from  “Ist-order  Navier- Stokes”  methods  to  DSMC.  These  alternative  solu¬ 
tions  attempt  to  extend  CFD  to  a  higher  order  approximate  solution  of  the  Boltzmann 
equation,  or  a  simplihed  version  thereof,  by  either  extending  the  computation  of  the 
dissipation  terms,  or  rewriting  them  in  terms  of  relaxation  terms,  and  treating  each 
mode  of  energy  separately. 

1.9  Bulk  Viscosity  and  Nonequilibrium 

It  is  useful  to  describe  what  we  mean  by  nonequilibrium  from  the  macroscopic 
sense.  Here  we  are  meaning  there  are  no  external  forces  on  the  gas,  and  thereby  no 
gradients.  In  the  absence  of  gradients  and  external  forces,  the  gas  is  in  an  overall 
inertial  rest  state.  No  spatial  changes  occur  (convective),  nor  are  there  any  changes 
with  advancing  time.  In  a  dilute,  monatomic  gas  this  is  fairly  simple.  In  a  polyatomic 
gas,  internal  energy  modes  exist,  which  are  excited  and  de-excite  between  quantized 
levels  through  collisional  processes.  That  is  to  say,  collisions  with  energy  comparable 
to  the  difference  between  quantized  energy  states  have  a  potential  to  trade  some 
translational  energy  for  the  internal.  These  collisions  are  called  “inelastic.” 

For  a  monatomic  gas,  the  Boltzmann  equation  (1.5)  can  be  used  as  a  starting 
point  for  the  microscopic.  At  rest,  there  should  be  no  rate  of  change  in  the  velocity 
distribution  function.  This  state  can  in  general  be  attained  in  a  gas,  when  the  right 
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hand  side  integrand  is  zero: 

r/r -//i  =  o  (1.7) 

This  relationship  is  one  of  collisional  invariance,  and  is,  in  fact,  closely  related  to 
the  formulation  of  the  equilibrium  distribution  function,  the  principle  of  detailed 
balance,  and  the  Boltzmann  H-theorem  (a  microscopic  version  of  the  second  law  of 
thermodynamics).  The  principle  of  detailed  balance  states  that  for  every  collision, 
there  is  a  corresponding  inverse  collision,  restoring  the  system  to  the  equilibrium 
distribution.  The  product  of  the  starred  quantities  is  known  as  the  “gain  term”, 
and  the  unstarred  the  “loss  term”.  It  is  interesting  to  note  that  quantum  mechanics 
had  not  been  developed  by  the  time  this  principle  was  first  derived.  The  principle  is 
somewhat  at  odds  with  the  Heisenberg  uncertainty  principle,  which  states  that  there 
is  a  finite  lower  bound  on  the  precision  one  can  observe  both  the  momentum  and 
position  of  a  particle. 

For  a  polyatomic  gas,  in  general  there  is  no  principle  of  detailed  balance.  The 
presence  of  the  internal  modes  complicates  matters  signihcantly.  Generalized  Boltz¬ 
mann  equations  have  been  formulated  to  try  to  describe  the  momenta  and  positions 
of  these  particles,  in  both  classical  [23],  and  quantum  mechanical  versions  [2:145]. 

The  continuum  description  of  the  fluid  can  be  obtained  from  the  Boltzmann  (in 
the  case  of  a  dilute  monatomic  gas)  [15],  or  generalized  Boltzmann  equations  (in  the 
case  of  a  polyatomic  or  electronically  excited  monatomic  gas)  [2:221-266].  The  form  of 
the  equations,  but  also  more  importantly  the  transport  coefficients  can  be  calculated 
using  the  information  contained  in  the  intermolecular  potential.  For  a  monatomic  gas, 
this  is  spherically  symmetric.  For  the  polyatomic,  it  is  not.  For  translation-rotation 
energy  exchanges,  this  anisotropy  is  very  important.  See,  for  example  [2],  [24]  for 
greater  detail. 

This  work  follows  that  of  Carr  [11],  and  in  his  work,  he  suggests  new  temperature- 
dependent  models  for  transport  coefficients,  and  the  inclusion  of  bulk  viscosity  in  the 
Navier-Stokes  Equations  to  extend  their  validity  in  altitude  in  the  continuum  and 
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near-continuum  regimes.  It  is  the  goal  of  this  work  to  treat  rotational  nonequilibrium 
in  the  Navier-Stokes  equations  as  a  way  of  improving  the  high- altitude  accuracy  of 
CFD  solutions.  This  work  does  not  attempt  to  solve  all  the  problems  of  both  transla¬ 
tional  and  rotational  nonequilibrium,  but  rather  to  gauge  the  effect  of  including  bulk 
viscosity  in  the  continuum  CFD. 

1.10  What  is  bulk  viscosity? 

Bulk  viscosity  is  a  resistance  to  a  volumetric  change  in  a  gas,  such  as  an  ex¬ 
pansion  or  contraction.  Expansions  and  contractions  are  mathematically  represented 
as  the  divergence  of  the  velocity  held.  Bulk  viscosity  linearly  relates  the  divergence 
to  a  difference  in  two  pressures:  thermodynamic  equilibrium  pressure  and  mechanical 
pressure,  or  one  third  of  the  trace  of  the  stress  tensor.  In  this  way,  any  nonequilibrium 
phenomena  that  makes  these  two  values  different  (rotation,  vibration,  other  internal 
mode  excitation)  could  be  treated  with  bulk  viscosity,  assuming  “small  deviations” 
from  equilibrium.  The  present  work  will  consider  only  the  rotational  contribution 
to  the  bulk  viscosity.  Traditional  Navier-Stokes  equations,  and  the  CFD  built  on 
them,  assume  that  the  rotational  modes  of  energy  are  always  in  equilibrium  with  the 
translational  (and  completely  neglects  other  internal  energy  modes).  Several  early 
works  state  that  the  equilibration  is  always  within  about  4-5  collisions,  and  therefore 
this  assumption  should  not  cause  a  problem,  with  the  exception  of  the  calculation  of 
shock  structure.  The  experimental  data  of  Alsmeyer  [25]  tends  to  agree  with  this. 
Other  experiments,  e.g.  Carnevale  et  al  [26],  Prangsma,  et  al  [27], Belikov  et  al  [28] 
and  some  models  [29],  [30]  disagree  with  this,  however.  Alsmeyer’s  data  was  gath¬ 
ered  from  density  prohles  in  Nitrogen,  at  room  temperature,  but  very  low  pressure 
(and  therefore  density).  Carnevale  et  al  used  ultrasonic  measurements  to  calculate 
rotational  relaxation  time,  and  Belikov  et  al  used  a  jet  expanding  into  near-vacuum 
conditions. 

Bulk  viscosity  can  be  found  in  the  Navier-Stokes  stress  term,  to  be  described 
in  detail  in  Chapter  II.  Vincenti  and  Kruger  [31:407-412]  and  McCourt  et  al  [2] 
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(as  volume  viscosity,  r]y,  as  is  common  in  the  Molecular  Physics  community)  give 
arguments  for  its  use.  It  should  be  noted  that  the  bulk  viscosity  treatment  has  been 
the  subject  of  controversy  for  quite  some  time  [32],  A  historical  treatment  of  bulk 
viscosity  (at  least  through  1999,  from  the  aeronautical  engineering  perspective)  can 
be  found  in  Argrow  and  Graves  [33].  Arguments  for  and  against  the  validity  of 
bulk  viscosity  as  a  physical  relaxation  parameter  can  be  found  in  Meador  [34],  and 
Emanuel  [35].  Some  of  the  controversy  is  over  the  apparent  frequency-dependence  of 
sound  attenuation.  Volume  viscosity,  as  it  is  seen  in  the  molecular  physics  literature  [2: 
280],  equation  (6.2-24),  is  related  to  bulk  viscosity,  as 
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where  on  is  the  sound  frequency,  rj  is  shear  viscosity,  a  is  the  absorption  coefficient. 
Cad  is  the  adiabatic  sound  speed  in  the  limit  of  zero  frequency,  A  is  the  thermal  con¬ 
ductivity,  k  is  Boltzmann’s  constant,  and  rjy  is  bulk  viscosity.  Using  sound  absorption 
experiments  makes  the  value  of  bulk  viscosity  calculated  from  these  kinds  of  exper¬ 
iments  very  sensitive  to  small  errors  in  the  measurement  of  a.  Also,  in  [2:284],  it  is 
explained  that  ultrasonic  sound  absorption  experiments  only  measure  the  rotational 
relaxation  part  of  bulk  viscosity.  The  reason  for  this  is  thought  to  be  the  much 
longer  relaxation  time  of  vibrational  modes,  and  the  high  frequency  of  the  sound. 
The  requirement  on  ultrasonic  measurements  of  absorption  coefficient  is  ±1%,  to  be 
considered  useful  for  measuring  bulk  viscosity.  Expansion  measurements  have  been 
subject  to  this  scrutiny  as  well. 


1.11  Scope  of  the  Thesis 

In  this  work.  Chapter  II  will  outline  the  history,  theory,  and  experimental  results 
for  rotational  relaxation,  the  basic  treatments  for  rotational  relaxation  within  the 
hydrodynamics  and  DSMC  frameworks.  The  recent  rotational  relaxation  calculations 
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of  Cappelletti,  Vecchiocattivi,  Pirani,  Heck,  and  Dickinson  [36:491]  are  discussed,  and 
used  as  the  basis  for  bulk  viscosity  in  an  axisymmetric  CFD  code,  AFIT-2D. 

Chapter  III  outlines  the  numerical  methods  within  AFIT-2D,  and  the  meth¬ 
ods  behind  a  one-dimensional  Navier-Stokes  solver  for  shock  structure,  based  on  the 
method  of  Gilbarg  and  Paolucci  [37],  used  to  show  the  effect  of  bulk  viscosity  on 
shock  structure.  It  should  be  pointed  out  that  translational  nonequilibrium,  be¬ 
yond  what  limited  amount  is  afforded  by  the  Navier-Stokes  equations,  is  beyond  the 
present  scope.  Use  and  modihcation  for  the  present  work  of  an  algebraic  grid  adaption 
routine,  borrowed  from  Gnoffo,  Hartung,  and  Greendyke  [38:2-3]  is  presented.  De¬ 
scription  of  the  two-dimensional  axisymmetric  test  cases  from  a  2003  code  validation 
experiment  [1]  is  described. 

Chapter  IV  presents  the  results  of  the  continuum  calculations,  comparing  the 
zero  bulk  viscosity,  constant  bulk  viscosity,  and  the  bulk  viscosity  calculated  from  the 
data  of  Cappelletti  et  ah  are  compared  in  both  ID  and  2D  cases.  Where  applicable, 
DSMC  data  from  Carr  [11]  is  used  for  comparison. 

Chapter  V  presents  conclusions  on  the  effect  of  including  the  bulk  viscosity,  and 
its  usefulness  in  turning  calorically  perfect  gas-based  CFD  into  a  preliminary  design 
tool  for  high-altitude  CFD.  Recommendations  for  future  work  are  presented. 
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II.  Rotational  Relaxation 


IN  classical  kinetic  theory,  rotational  energy  is  expressed  in  terms  of  ‘degrees  of 
freedom.’  Here  the  number  of  axes  about  which  the  molecule  can  rotate,  and  posess 
a  signihcant  amount  of  energy  is  multiplied  hj  Ths/^,  which  amounts  to  a  maximum 
entropy  state.  This  statement  of  course,  is  made  assuming  that  the  numbers  of  states 
available  is  continuous  and  proportional  to  temperature. 

The  typical  continuum  stress  tensor  for  a  Newtonian  fluid  is  given  by 


P^ij  T  p 
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where  A  is  commonly  labeled  the  second  coefficient  of  viscosity,  which  is  related  to 
the  bulk  viscosity,  ps  by  the  dehnition 


A  =  Ps  —  -p 


(2.2) 


Typically,  a  “mechanical  pressure”  is  dehned  as  the  average  of  the  dilatational  stresses  [2: 
123],  [31:410]: 

1  duh 

P  =  =  (2.3) 

The  “thermodynamic  pressure”  is  dehned  by  the  equilibrium  pressure  for  the  given 
temperature: 

p  =  pRT  =  nksT  (2.4) 

where  n  is  the  number  density,  ks  is  Boltzmann’s  constant.  The  two  representations 
are  related  by  p  =  mn  and  R  =  kB/m,  where  m  is  the  molecular  weight.  From  kinetic 
theory  (see  Vincenti  and  Kruger  [31:410]  or  McCourt  et  al  [2:122-123]): 

1  —  c)ii  ■ 

P-P  =  P-^PC^  =  PB-^^  (2.5) 

therefore  if  p  =  p,  ^  ps  =  0  Thermal  velocity,  Ci,  is  dehned  as  the  molecular  speed 
Ci  minus  the  mean  how  velocity  Ui. 
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It  should  be  mentioned  that  taking  internal  energy  modes  into  account  in  the 
energy  equation  was  hrst  treated  by  Eucken  [2:62,124],  This  expression  for  Eucken 
factor  relates  the  shear  viscosity,  isochoric  (or  constant  volume)  heat  capacity,  and 
thermal  conductivity: 


Ie  — 


k 


km 

fJjCy 


(2.6) 


This  result  should  look  familiar  to  aerodynamicists  as  [cp/cv)/ Pr.  For  monatomic 
gases,  this  value  is  5/2.  For  polyatomic  gases,  approximate  empirical  relationships 
are  often  used  [2:124]: 

(2.7) 


t 


fE  =  fE-+f, 


c 

Cl  V 


where  i  denotes  internal  energy,  and  t  denotes  translational  energy,  and  [39:30]: 


Pr 


47 

7.O87  -  1.80 


(2.8) 


Stokes’  Hypothesis  sets  the  bulk  viscosity  to  zero,  which  effectively  sets  the 
two  pressures  always  equal.  This  hypothesis  therefore  implies  equilibrium  condi¬ 
tions  throughout  the  flow.  It  is  experimentally  and  analytically  verihed  for  dilute 
monatomic  gases.  By  dilute  here  we  also  are  bringing  in  the  idea  of  a  low-density 
limit  of  the  bulk  viscosity.  Even  for  a  monatomic  gas,  the  bulk  viscosity  is  nonzero 
when  the  density  is  not  sufficiently  low.  See  Rah  and  Eu  [40]  for  more  detail.  In 
fact,  the  other  transport  properties  used  in  the  Navier-Stokes  description  in  terms  of 
temperature  alone  are  the  low-density  limit  transport  properties,  and,  in  general,  are 
density  dependent  [39:25-31]. 

The  higher  density  gases  also  cease  to  be  treatable  by  the  ideal  gas  law 


p  =  pRT  =  nksT  (2.9) 

but  are,  rather  by  the  van  der  Waals  equation  of  state  [41:22] 

(p+^){v-b)  =  RT  (2.10) 
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which  is  also  known  as  the  real  gas  equation  of  state.  Bulk  viscosity  has  its  greatest 
effect  in  non-equilibrium  situations,  especially  places  like  shocks  where  the  transla¬ 
tional  modes  of  energy  are  the  first  to  be  excited,  and  there  is  a  collisional  process  of 
exciting  internal  energy  modes.  In  a  monatomic  gas,  kinetic  theory  predicts  zero  bulk 
viscosity,  as  until  quite  high  temperatures,  there  are  no  internal  modes  to  be  excited, 
and  the  only  kind  of  nonequilibrium  present  in  this  situation  is  translational. 

In  a  diatomic  gas,  such  as  N2  or  air,  the  rotational  mode  can  lag  the  transla¬ 
tional,  but  only  by  a  few  collisions.  So,  in  a  dense  (say  sea-level  standard  conditions) 
environment,  the  region  of  non-equilibrium  in  a  supersonic  flowfield  is  very  small  and 
can  often  be  neglected.  However,  at  high  altitudes,  in  lower  density,  this  effect  can 
be  more  pronounced. 

The  idea  of  using  bulk  viscosity  as  a  way  of  compensating  for  small  degrees  of 
nonequilibrium  is  not  new  [15],  [31].  A  paper  by  Graves  and  Argrow  [33]  outlines 
a  history  of  the  bulk  viscosity,  including  its  use  to  predict  nonequilibrium  behavior. 
The  bulk  viscosity  has  been  reported  to  be  frequency  dependent  by  some  in  the  field 
of  acoustics.  However,  a  search  of  the  literature  reveals  little  quantitative  information 
about  what  a  “small  degree”  of  nonequilibrium  really  is.  Heck  et  al.  remark  that 
there  is  a  lack  of  data  on  the  bulk  viscosity,  and  that  obtaining  it  is  difficult. 

Values  of  bulk  viscosity  due  to  Heck  et  al.  [24]  [42]  will  be  used  in  this  study. 
These  are  based  on  collision  cross-sections  calculated  from  a  relatively  new  model  of 
the  intermolecular  potentials  for  nitrogen.  These  potentials,  unlike  the  hard  sphere. 
Maxwell  molecule,  Sutherland,  and  Lennard-Jones,  are  anisotropic,  since  Nitrogen 
molecules  are  not  spherically  symmetrical.  This  anisotropy  of  the  potential  determines 
how  rotational  energy  is  transferred  between  colliding  molecules. 

2.1  The  Boltzmann  Equation  for  Dilute  Monatomic  Gases 

The  Boltzmann  equation  can  be  derived  from  the  6N-degree  of  freedom  Liouville 
equation  for  the  probability  distribution  function  of  a  molecule’s  position  and  velocity 
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in  an  N-molecule  system.  This  treatment  is  beyond  the  scope  of  this  thesis,  and  can 
be  found  in  Sone  [43:257-269].  The  Liouville  equation  itself  is  a  result  from  the  held 
of  statistical  mechanics. 

It  is  somewhat  a  misnomer  to  call  a  physical  gas  dilute.  A  dilute  gas  is  a  math¬ 
ematical  construction,  coming  from  an  assumption  that  no  more  than  two  molecules 
interact  (collide)  at  the  same  time  [2:2].  This  construction  is  also  known  as  the  binary 
collision  model.  The  effect  of  other  molecules  during  “free  hight”  and  “collisions”  are 
ignored.  No  matter  the  form  of  the  intermolecular  potential,  that  potential  is  treated 
as  vanishing  beyond  some  hnite  distance.  In  the  physical,  it  is  practically  impossible 
to  verify  that  this  is  ever  the  case.  A  dilute  gas  then  could  be  dehned  as  one  for  which 
predicted  macroscopic  behavior  using  the  dilute  assumption  “agrees  well”  with  mea¬ 
sured  macroscopic  behavior.  That  is  to  say,  one  for  which  the  amount  of  non-binary 
collisions  has  a  negligible  effect  on  the  overall  compared  values. 

In  a  gas  at  rest  and  in  thermal  equilibrium  with  its  boundaries,  the  veloci¬ 
ties  of  the  molecules  tend  to  be  distributed  according  to  a  Maxwellian  Distribution 
distribution  [31:35]: 


f{a)  =  f{c,,c2,c,) 


(  m  \3/2 


m 


2kT 


(C^  +  C?_+Cl) 


(2.11) 


Recall  the  Boltzmann  Equation  from  chapter  I  [9:53]: 

—  (n/)  +  c-— (n/)  +  F-— (n/)  =  /  /  n^{f*  -  f  fi)CradVLdci,  (2.12) 

»/  —  oo  V  0 

where  /  is  the  pre-collision  velocity  distribution  function  of  the  molecule  selected  for 
collision,  /i  is  the  pre-collision  velocity  distribution  function  of  the  collision  partner, 
and  the  starred  versions  are  the  post-collision  representations  of  the  same,  c  is  the 
molecular  velocity  of  the  molecule  associated  with  /.  F  is  an  external  force  (such  as 
applied  electric  held). 
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Chapman  and  Enskog  separately  developed  an  expansion  about  the  normal 
distribution,  here  as  equation  2.11  [15:107]: 

/  =  /'“>  +  /«  +  /'^>  + 

Enskog  introduced  a  parameter  6,  which  can  be  related  to  Knudsen  number  [15:115]: 

/  =  +  •  •  •  (2.14) 

This  process  ends  in  an  expansion  of  unknown  functions  0"^  and  Knudsen  number: 

/  = /(°)(l  +  77n0i  +  i7n202  +  ...)  (2.15) 

Lumpkin  [21:15]  points  out  this  expansion,  used  in  the  Boltzmann  equation  can  be 
shown  to  give  the  Euler  equations  when  truncating  at  0  order,  the  Navier-Stokes 
equations  when  carried  to  first  order,  and  the  Burnett  equations  when  carried  to 
second  order,  but  in  general,  is  asymptotic  or  divergent.  Results  can  be  obtained 
for  small  Kn,  but  the  expansion  exhibits  nonconvergent  behavior  for  Kn  3>  1.  This 
behavior  is  often  cited  as  the  reason  why  the  continuum  assumption  fails  for  large 
Kn. 

2.2  Classical  Rotational  Relaxation  Results 

The  rotational  collision  number  Zr  is  a  measure  of  the  ratio  of  rotational  relax¬ 
ation  time  r,  the  time-constant  for  rotation  to  equilibrate  with  translation  in  a  heat 
bath,  to  mean  time  between  collisions  Tc'. 

Zr  =  -  (2.16) 

The  rotational  relaxation  time  is  essentially  a  time  constant  in  a  simple,  first- 
order  ordinary  differential  equation,  originally  proposed  by  Jeans  [44],  as  cited  by 


(2.13) 
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Lumpkin  [21:52]: 


(2.17) 


dCr  E^^(Tt)  —  Er 

dt  T 

This  equation  has  the  same  form  as  the  Landau- Teller  equation  for  vibrational  re¬ 
laxation.  It  assumes  translational  equilibrium,  meaning  the  velocity  distribution  of 
molecules  during  this  process  is  Maxwellian. 

Parker  [29:454]  in  his  classic,  well  known  article,  solved  a  system  of  classical 
equations  of  motion  analytically  to  obtain  an  expression  for  the  rotational  collision 
number,  Zr'. 

^  _ _ ^oo _ 

l  +  iri(S/r)J  +  (f +  2)  (S/T) 
where  Z^o  is  the  maximum,  or  perhaps  more  appropriately,  asymptotic  rotational 
collision  number,  and  S  is  known  as  Sutherland’s  Constant  by  White  [39:28].  It  is  a 
measure  of  the  “well  depth”  [29:458]  of  the  intermolecular  potential.  Note  that  the 
coefficient  of  the  {S/T)  term  has  the  correction  of  Brau  and  Jonkman  [30:481]  applied. 
One  assumption  of  this  work  was  that  of  a  gas  with  initially  nonrotating  molecules. 
Another  was  one  of  “impulsive”  or  “sudden”  collisions,  where  the  duration  of  an 
intermolecular  collision  was  treated  as  much  shorter  than  the  rotational  period  of  the 
molecule. 

This  form  of  rotational  collision  number  is  still  in  use  in  DSMC  today  [45:2]. 
Brau  and  Jonkman  re-solved  this  problem  stochastically,  hrst  by  removing  Parker’s 
hrst  assumption,  and  obtained  exactly  one-half  his  result  [30:481].  Then,  they  came 
to  a  solution  assuming  “adiabatic”  collisions,  that  is  to  say,  collisions  for  which  the 
rotational  period  is  much  shorter  than  the  collision  time,  and  found  that  for  this  case, 
the  idea  of  a  single  relaxation  time  is  no  longer  valid. 

2.3  Phenomenological  and  Semiclassical  Results 

Some  insight  into  the  dehnition  of  phenomenological  model  can  be  obtained 
from  Bird  [46:2,6].  The  premise  is  a  model  which  qualitatively  agrees  with  phenomena 
observed  in  nature.  The  opposite  end  of  the  spectrum  is  the  idea  of  ab  initio  or  “from 


(2.18) 


22 


first  principles.”  In  the  case  of  molecular  physics,  the  principles  of  quantum  mechanics 
are  treated  as  the  first  principles.  In  classical  treatments,  Newton’s  Laws  are  the  hrst 
principles. 

If  the  quantum  mechanics  are  the  more  appropriate  starting  point  in  fluid  flow, 
then  the  Navier-Stokes  (here  the  meaning  is  the  momentum  conservation  equations 
for  which  Navier  and  Stokes  are  strictly  responsible)  constitutive  models  are  a  very 
successful  example  of  phenomenological  modeling.  Quantum  mechanics  were  not 
available  at  the  time,  and  microscopic,  kinetic  theory  models  were  not  used  in  their 
derivation. 

The  DSMC  codes  produced  by  Bird  [9]  use  the  Larsen-Borgnakke  phenomeno¬ 
logical  model  for  internal  energy  relaxation.  This  model  uses  a  constant  rotational 
collision  number,  which  for  Nitrogen  has  a  default  value  Zr  =  5  [47].  The  DSMC 
code  available  to  Carr  [11]  was  MONACO,  from  Boyd’s  group,  which  uses  Parker’s 
model  [45:2]  for  rotational  relaxation.  This  model  uses  a  rotational  collision  number 
of  the  form  in  equation  (2.18).  The  default  values  of  S  and  Z^q  are  those  of  the 
original  Parker  model.  Carr  [11:41]  uses  Z^o  =  15.7,  S  =  80K  in  his  solutions.  The 
rotational  collision  number  is  related  to  the  variable  hard  sphere  model,  as  outlined 
in  Lumpkin,  Haas,  and  Boyd  [48].  Wysong  and  Wadsworth  give  an  assessment  of 
phenomenological  rotational  relaxation  models  for  DSMC  [49]. 

2.4  Generalized  Boltzmann  Equation  Rotational  Relaxation  Results 

Cappelletti  et  al  [36]  compute  a  potential  energy  surface  PES,  modifying  existing 
work  by  Heck  and  Dickinson  [24] ,  [42] ,  and  then  compute  transport  coefficients  using 
Heck  and  Dickinson’s  semiclassical  solver  based  on  a  second-order  Chapman-Enskog 
expansion  of  the  generalized  Boltzmann  equation  of  Curtiss,  Kagan,  and  Maksimov 
as  described  in  [2:145-158].  One  interesting  feature  of  this  work  is  a  comparison 
between  the  new  PES  and  an  earlier  surface,  the  AWJ  potential.  The  AWJ  potential 
was  obtained  from  an  ab  initio  short-range  interaction.  The  new  PES,  deemed  PES8 
was  a  modihcation  of  the  earlier  AWJ  potential,  adjusted  to  better  match  observed 
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relaxation  cross-sections.  Both  PES  are  of  the  form 


l/(R,ri,r2)  =  (47r)^/^  ^  ri,  ra)  (2.19) 

La,Lb,L 

where  A  is  a  set  of  orthormal  functions,  R  is  the  vector  from  the  center  of  molecule 
A  to  that  of  molecule  B,  and  ri  and  r2  [36:486]. 

The  article  of  Cappelletti  et  al  [36]  does  not  explicitly  discuss  the  form  of  the  hy¬ 
drodynamics  equations  which  come  from  this  particular  application  of  the  Chapman- 
Enskog  expansion.  The  present  work  will  use  the  results  for  fis  in  the  Navier  Stokes 
equations.  The  incorporation  of  bulk  viscosity  in  the  Navier-Stokes  equations  in  this 
way  could  be  considered  “phenomenological  modeling”  by  the  dehnition  of  Bird  [9]. 

It  should  be  noted  that  the  aerodynamics  community  is  just  starting  to  take 
interest  in  generalized  Boltzmann  equations,  for  example,  see  Cheremissin  and  Agar- 
wal  [50].  These  workers  report,  as  does  [2],  the  Wang-Chang-Uhlenbeck-de  Boer  [51] 
equation,  seen  in  the  aerodynamics  literature  in  the  1950s  and  1960s,  does  not  han¬ 
dle  properly  the  degeneracy  of  rotational  states.  Cheremisin  and  Agarwal  present 
one  dimensional  shock  solutions,  taking  on  the  order  of  weeks  to  compute  using  the 
generalized  Boltzmann  equation  they  have  been  using. 


2.5  Bulk  viscosity  from 

The  bulk  viscosity  is  related  to  the  the  rotational  relaxation  time,  r,  as 


-pr 


(2.20) 


where  c„  is  the  isochoric  or  constant-volume  heat  capacity,  and  Crot  is  fhe  isochoric 
rotational  heat  capacity  [36].  The  shear  viscosity  coefficient  is  dehned  by  Parker  [29] 
to  be  related  to  the  mean  time  between  collisions,  Tc,  by 


Tc  = 


Tip, 

4p 


(2.21) 
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The  rotational  collision  number,  is  defined  by  equation  2.16.  Solving  equations 
(2.16),  (2.20)  and  (2.21)  for  /isZ/x,  we  obtain 


TT 

4 


(2.22) 


In  the  present  work,  Crot  =  ks  is  assumed.  It  should  be  mentioned,  however,  that 
correlations  for  Crot  are  available  [36:491].  This  assumption  allows  a  simpler  expression 
for  the  ratio  of  bulk  to  shear  viscosity: 


=  ^(7  - 


(2.23) 


Zr  is  fit  to  the  reported  values  from  [36:492]  with  a  polynomial  in  temperature: 


Zr{T)  =  fsT^  +  +  hT  +  /o  (2.24) 


where  /s  ~  1.034E  —  9,  /2  ~  — 3.835E  —  6,  /i  ~  1.104E  —  2  and  /o  ~  1.905.  This 
fit,  along  with  the  reported  data,  is  depicted  in  Figure  2,  along  with  the  models  due 
to  Parker  [29]  and  [30]  for  comparison.  These  models  are  commonly  used  in  DSMC. 
Here,  the  Parker  model  has  its  original  values  of  Z^o  and  S.  Note,  the  newer  work 
reports  a  higher  rotational  collision  number,  and  therefore  higher  bulk  viscosity.  Also, 
though  the  data  from  Cappelletti  et  al  is  not  reported  above  about  1400K,  it  does 
not  suggest  an  asymptotic  value  at  higher  temperatures.  The  collision  model  used  for 
determining  Z^  uses  the  infinite-order  sudden  approximation.  This  approximation  is 
not  valid  at  high  collision  energies,  and  therefore  is  limited  to  this  low  temperature 
range.  For  more  detail  on  the  difference  between  sudden,  adiabatic,  and  close-coupled 
approximations,  the  reader  is  referred  to  [2,49,52]. 
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Translational  Temperature,  K 


Figure  2:  Rotational  collision  number 
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III.  Numerical  Methods 


The  basic  framework  for  testing  the  bulk  viscosity  in  the  continuum  model  is 
AFIT-2D,  developed  at  the  Air  Force  Institute  of  Technology.  This  code  being 
small,  flexible,  and  with  all  source  code  available,  was  the  natural  choice  for  this  work. 


3.1  History  of  AFIT-2D 

AFIT-2D  is  a  hnite  volume  Navier-Stokes  code  for  simulation  of  two  dimen¬ 
sional  compressible  flows.  The  present  work  uses  an  augmented  version  for  axisym- 
metric  flow  that  retains  the  bulk  viscosity  in  the  viscous  terms.  The  code  has  several 
user-selectable  time  integration  schemes  and  flux  schemes.  The  flow  variables  and 
transport  coefficients  are  nondimensionalized  by  Mach  number,  Reynolds  number  and 
Prandtl  number.  Recall  from  Chapter  II  that  the  use  of  a  constant  Prandtl  number 
is  equivalent  to  Eucken’s  approximation  for  heat  transfer.  While  the  solver  is  written 
in  Reynolds  Averaged  Navier  Stokes  form,  which  can  compute  eddy- viscosity-based 
turbulent  flows,  all  solutions  presented  in  this  work  were  computed  as  laminar.  The 
Reynolds  numbers  of  the  flows  considered  were  all  well  within  the  laminar  flow  regime, 
as  conhrmed  by  experiment  [1]. 


3.2  Methodology  of  AFIT-2D 

Axisymmetric  flow  is  included  as  a  special  case  by  means  of  a  source  term  in  the 
Navier-Stokes  equations.  The  formulation  in  conserved  variable  vector  form  adapted 
from  Hoffman  and  Chiang  [53:Vol  H,  456-457]  is 


dQ 
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dE 
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where  a  =  0  for  the  case  of  a  2D  cartesian  problem,  a  =  1  for  an  axisymmetric 
problem, 
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where  A  is  defined  by  equation  (2.2). 

In  typical,  calorically  perfect  gas-based  CFD, 
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Sutherland’s  model  for  shear  viscosity  is  often  used  in  CFD: 

(3.12) 

f^ref  \^ref  /  1  -r  O 

The  grid  transformation  presented  by  Gnoffo,  Hartung  and  Greendyke  [38]  will 
be  used  to  cluster  cells  in  the  vicinity  of  the  shock  and  the  wall,  where  gradients  are 
highest.  Bulk  viscosity  will  be  first  calculated  from  curve-fit  equation  2.24. 


3.3  Calculation  of  Convective  Flux 


For  convective  flux  calculation  in  AFIT-2D,  the  flux  difference  splitting  method 
of  Toro,  Spruce  and  Speares  [54]  (HLLC)  is  used,  selecting  wavespeeds  as  described  in 
the  paper  by  Batten  et  al  [55].  This  scheme  restores  the  contact  discontinuity  to  the 
HLL  scheme  of  Harten,  Lax,  and  van  Leer  [56].  The  HLLC  flux-difference  splitting 
method  is  defined  in  [55]  as 
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The  recommended  wave  speeds  from 


Fl 

:  ^L>0 

Ft  : 

:  <  0  <  A* 

Ft 

:  S*  <{)<  Sr 

Fr  : 

'■  Sr  <0 

[55]  are  used: 


(3.13) 


Si  =  mm[A,(l7,),A,(l7'^); 


(3.14) 
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and 

Sr  =  max  [XmiU^n,  ^m{Ur)]  (3.15) 

Also, 

=  Q  (3.16) 

\i  =  Ui  —  ai  ^  \i{Q)  =  Urob  —  OjRob  (3.17) 

and 

\m  =  U  +  a  (3.18) 

Starred  values  are  defined  as 


F:  =  F,  +  S,{Q*-Qi),  z  =  L,R 


(3.19) 


where  Fr,  Fr  are  the  exact  fiuxes  in  cells  adjacent  to  the  face  of  interest.  Rearranged, 
this  is  the  same  as 

S,Q:  -F:  =  S,Q,  -  F,  (3.20) 

which  for  the  first  term  in  Q  and  F,  was  demonststrated  as 


-  p:ui^r=l) 

pKs,-s*) 

Pi 


Si  Pi  PiUi 

P^{Si  -  ij^) 
Si-iii 


(3.21) 


The  other  lines  are  substituted  into  (3.20)  as  well.  For  M-momentum: 


^  77*  =  S* 

Si{piUi)*  -  {piUiYpf-  p*n^ 


=  Si{piUi)  -  {piUi)Ui  -  PiTl^ 


{PiUi}*  = 


Ui  = 


[PiUi 


Ui  ,  ip*  -  Pi)n:, 


Si  -  S* 


+ 


R  - 


{piUi){Si  -  Ui)  {p*  -pi)n^ 


pns,-s*)  pm-s*y 


(3.22) 

(3.23) 
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Substituting  (3.21)  into  (3.23), 


U*  =  Ui  + 


{p*  -Pi)n^ 

P^{S,  -  U) 


Similarly, 


{piViY  =  {piVi 


Si-Ui  ^  {p*  -  Pi)ny 


S,-S 


&  - 


or 


Now,  for  the  fourth  line: 


jp*  -  Pi)ny 
*  pYSi  -  U)  ■ 


S,El  -  {El  +p*)ir  =  S,E,^  -  {E,^  +pYU, 


El{Si-  S*)  -P*S* 

e; 


EtYs^-uY-p^u, 

s,-u,  p*s*-p,u, 
Si-S*^  s,-  s* 


(3.24) 


(3.25) 


(3.26) 


(3.27) 


S*  is  dehned  by 

_  jj*  _  PrUr{Sr  —  Ur)  —  PlUl{Sl  —  Ul)  +  Pl  —  Pr 
““  Pr{Sr-Ur)-pl{Sl-Ul) 


p*  is  given  by 

p*  =  pYu^-sY{u,-s*)+p^. 


(3.28) 


(3.29) 
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The  Roe- Averaged  Variables  are  defined  by 


where 


R 

- 

(3.30) 

V  PL 

PRoe 

=  RpL 

(3.31) 

^Roe 

_  Rur  +  Ul 
~  I-PR 

(3.32) 

VRoe 

_  Rvr  +  Vl 

1  +  R 

(3.33) 

Uros 

^Roe^x  “1“  ^Roe^y 

(3.34) 

^ORoe 

_  Rho^  +  ho^ 
l  +  R 

(3.35) 

2 

^Roe 

=  (7-1)  Voe  -  ^(“Le  +4oe) 

(3.36) 

Et+p 

P 


(3.37) 


3.4  Calculation  of  Viscous  Flux 

3.4-1  Gradients  -  Green-Gauss  Approach.  The  procedure  for  calculating  a 
flux  involves  a  diamond-shaped  control  volume,  shown  in  Fig.  3(for  an  i-face).  while 
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Figure  3:  Control  volume  for  viscous  flux  at  an  face  of  constant  i 


that  for  a  j-face  is  depicted  in  Fig.  4.  For  i-faces,  1  and  3  are  cell  values,  and  2  and 
4  are  node  values.  For  j-faces,  1  and  3  are  node  values,  and  the  2  and  4  faces  are  cell 
values.  Repeated  here  is  the  Green-Gauss  Theorem: 


f  dF2  dF^ 
Js  dx  dy 


Fidx  -|-  F2dy 


(3.38) 
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X2 


Figure  4:  Control  volume  for  viscous  flux  at  an  face  of  constant  j 

With  Fi  =  0,  F2  =  G  (G  is  a  property  within  the  flow  such  as  u,  v  or  T)  and 
approximating  the  right  hand  side  with  constant  values  of  u  at  each  face,  and,  on 
realizing  that  in  2  dimensions  S  =  V_,  (3.38)  becomes 


dx—^  ^  '-^/ace,aDg 
faces 


^  ^  Gi 


(3.39) 


where  V_  =  4(Aa;i3A|/24  —  Aa;24A|/i3)  Evaluating  on  the  dotted  control  volume,  this 
expression  yields 


dx 


1 

V 


-(Gi  +  G2)(|/2  —  Pi)  +  2(^2  +  G3)(|/3  —  1/2)  +  2^^^  ^4)(|/4  —  Ps)  +  2^^^  ~  ^t) 

(3.40) 


Expanding,  obtain: 


+  G1P2  —  G2P1  +  Q'fPfl  —  +  G2I/3  —  G3I/2  +  +  G3I/4  —  G4I/3 

~  +  G4I/1  +  —  G1I/4]  (3-41) 


Simplifying, 


(9a; 


Gi)A|/24  ~  (G4  —  G2)A|/i3 


(3.42) 
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With  Fi  =  G,  F2  =  0  and  approximating  the  right  hand  side  with  constant  values  of 
u  at  each  face,  and  realizing  that  in  2  dimensions  S  =  V_,  (3.38)  becomes 


dGi  dG  . 

dx  ^  j  ^ face face 

faces 


and,  by  the  same  process  as  above. 


dG  1 


”(^3  ~  Gi)Aa;24  +  (G4  —  G2)Aa;i3] 


dy  21^ 

On  substituting  flow  variables  for  (3.44),  (3.42), 

du  1  r. 


and 


dx  2V 

du 

1  r 

dy  ^ 

2V  ^ 

dv  1 

dx 

20 

dv 

1  r 

dy  ^ 

20  ^ 

dT  1 

dx 

^  20 

dT 

1  r 

dy  ^ 

20  ^ 

/13 


/13 


(3.43) 

(3.44) 

(3.45) 

(3.46) 

(3.47) 

(3.48) 

(3.49) 

(3.50) 


The  reader  is  reminded  /xs  =  A  —  2/3/r.  Also,  this  expression  is  divided  by 
reference  Reynolds  number,  and  the  substitution 


k 


1 

(7-1) 


k-L  kT  \ 

Ptl  Ptt  ) 


(3.51) 
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is  made,  where  for  the  present  work  ht  =  ^ 

Now,  if  G  represents  a  property  within  the  domain,  such  as  T,u,v,n,  then  we 
can  hnd  a  relationship  for  node  values  in  terms  of  an  average  of  surrounding  cells. 
Figure  5  contains  a  diagram  for  reference. 


{i  -  l,j  +  1) 


ihj  +  1) 

I 


(i  +  lj  +  l) 


{i  - 


x(f  -  1,  j) 

A  M 

k  A 

w 

-  l,i  -  1) 

w 

{hi  -  i)x 

(i  +  l,j) 


t 

ihJ  -  1) 


Figure  5:  Viscous  flux  stencil  diagram 


One  way  to  get  the  value  at  a  node,  which  will  be  used,  is  by  simple  averaging 


Gnode,i,j  ^i,G cell,i—l,j  G "b  Gcell,i—l,j  —  l  Gcell,i,j  —  l)  (3.52) 


Now  for  an  i-face,  in  this  case  the  left  face  of  cell  (i,  j),  by  Fig  3: 


Gi 

(3.53) 

G2 

G  node^ij 

(3.54) 

Gs 

Gcell^iJ 

(3.55) 

G4 

Gfiode^iJ+l 

(3.56) 

Rlld  Gi—jace^iJ 

=  x(G2  +  G4) 

(3.57) 
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For  a  j-face,  in  this  case  the  bottom  face  of  cell  {i,j) 


Cri  Gnode,i,j 
G2  =  Gcell,i,j-l 


G3 


G 


node,i+l,j 


G4 


and  G 


j-face,i,j 


-(Gi  +  G3) 


(3.58) 

(3.59) 

(3.60) 

(3.61) 

(3.62) 


Now,  we  can  write  the  viscous  fluxes  as  dehned  in  §2.4  of  Blazek  [57],  trading  Ux  and 
Uy  for  Ax  and  Ay,  since  area-scaled  fluxes  are  used  in  AFIT-2D; 


0 

AxTxx  "F  AyTxy 
-^x^xy  4“  ■^y'^yy 
^x^x  +  ^y^y 


(3.63) 


All  information  needed  for  either  an  i-face  or  j-face  in  (3.45)-(3.62)  is  now  available. 


3.5  Calculation  of  Time  Step 

The  time  step  of  the  simulation  in  the  case  of  an  explicit  time  integration  scheme 
such  as  the  second  order  Runge-Kutta  total  variation  diminishing  scheme  is  computed 
in  cell  {i,i)  according  to 


= 


+  A^)jj  +G(kl  +  kl)i,j 


(3.64) 


where  jg  the  volume  cell  (f,j).  A^  is  the  maximum  eigenvalue  of  the  convective 
flux  Jacobian  (the  spectral  radius). 
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3.6  Viscosity  in  the  Implicit  Operator  of  the  LUSGS  Scheme 

The  implicit  time  integration  scheme  of  Yoon  and  Jameson  [58]  is  employed  for 
this  work.  In  this  scheme  the  implicit  operator  is  approximately  factored  in  three 
parts: 

{D  +  L)D~\D +  U)AW^  = -R^  (3.65) 

where  D  represents  the  diagonal  terms  of  the  operator,  L  represents  the  lower  triangu¬ 
lar  terms,  U  represents  the  upper  triangular  terms,  AW'^  represents  the  change  in  the 
solution  at  time  level  n,  and  represents  the  residual,  or  explicit  spatial  operator. 

For  the  Euler  equations,  the  D~^  term  is  computed  as: 


D-^ 


1 


M  +  ^c) 


(3.66) 


(the  identity  matrix  is  implied/applied  in  the  code). The  average  areas  are  used  here 
and  below  for  the  visous  case,  in  D~^.  The  viscous  form  is  that  of  equation  (6.53)  in 
Blazek  [57:206].  First,  computing  viscous  spectral  radii: 


A.  = 


At  = 


/  4  7 
max  — ,  — 

Vsp’p 


max 


4  7 
3p’  p 


hL  PT 


Pr-L 

Pl 


+ 


PttJ  V 

(A^f 


Pt 


Pr  T,  Pr 


(3.67) 

(3.68) 


where  the  code  has  taken  7  >  4/3  as  hard-coded.  When  adding  bulk  viscosity,  4/i/3 
immediately  becomes  4/i/3  -|-  ps-  From  this  viscous  spectral  radius,  D~^  becomes: 


DA  = 


A  lam 


1  +  2DPJAIA'  +  AiAJ) 


(3.69) 


By  equations  (6.51)  in  Blazek  [57:205],  and  the  spectral  radius  simplihcation,  LAQ 
and  UAQ  become 


L  AC^  LAQiam  T  A^AQi—\  j  -\-  A^  AQi  j—\ 


(3.70) 
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and 

Ul^Q^isc  =  LAQiam  -  K^Qi+ij  -  AiAQi,,+i.  (3.71) 

3.7  Grid  Adaption 

To  align  the  cells  with  the  shock  the  procedure  of  Gnoffo,  Hartung,  and  Greendyke 
is  used  [38:2-3].  This  simple,  structured-grid,  algebraic,  line-by-line  procedure  as¬ 
sumes  one  of  the  curvilinear  grid  coordinates  moves  from  the  body  surface  to  the 
freestream.  One  transformation  is  used  to  resolve  the  boundary  layer,  based  on  a 
user-input  Reynolds  number  based  on  cell  height.  A  second  transformation  is  per¬ 
formed  to  normalize  the  far  held  “away  from  body”  distance  equal  to  one.  The  shock 
location  is  determined  by  sensing  a  user-input  fraction  above  the  freestream  density, 
pressure  or  temperature,  marching  from  the  freestream  to  the  body.  The  hnal  trans¬ 
formation  is  made  to  pull  grid  points  to  this  shock  location.  The  free  stream-side 
outer  boundary  is  placed  a  user-specihed  distance  beyond  the  sensed  shock  location. 
This  procedure  is  carried  out  for  every  value  of  non- surface-normal  grid  coordinate, 
on  the  face  centers,  and  the  nodes  are  placed  by  simple  averages  of  those  face  center 
positions. 

To  modify  this  procedure  to  accommodate  concave  geometries  (such  as  will  be 
used  in  this  work,  see  section  3.10,  rather  than  average  the  face  centers  directly,  an 
additional  constraint  restricts  motion  of  the  nodes  to  the  line  on  which  they  were  orig¬ 
inally  placed.  This  option  is  specihed  in  the  input  hie  as  “preserve  vectors  :  :  1.” 
With  this  option,  initial  grids  were  constructed  so  as  to  avoid  constant-r]  node  lines 
crossing  and  causing  a  grid  “folded  over”  on  itself.  This  subroutine,  in  Fortran  95, 
can  be  found  in  Appendix  A. 

3.8  Boundary  Conditions 

The  boundary  conditions  typically  used  for  Navier-Stokes  GFD  involve  no  slip  at 
the  wall,  and  no  pressure  gradient  at  the  wall.  Then,  the  only  property  left  to  specify 
is  the  temperature  gradient,  or  temperature  of  the  wall.  For  the  blow-down  tunnel 
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experimental  runs,  the  constant  wall  temperature  boundary  condition  is  appropriate. 


dp 

dn 


0 

0 


=  const 


(3.72) 

(3.73) 

(3.74) 


In  a  cell-centered  finite  volume  code  such  as  AFIT-2D,  the  boundary  conditions 
are  not  specified  directly  at  the  surface,  but  are  rather  specified  at  a  “ghost  cell” 
outside  the  flow  domain.  The  cell  interface  is  the  wall,  and  the  ghost  cell  properties 
can  be  set  such  that  the  flux  scheme  will  match  the  condition  of  no  mass  flow  through 
the  cell  face,  or  in  this  case,  calculate  zero  flow  velocity  at  the  cell  face.  In  AFIT-2D, 
the  nondimensional  pressure,  temperature,  and  temperature  are  related  by 

T  =  —  (3.75) 

P 

To  set  the  pressure  gradient  to  zero  is  simply  done  by  copying  the  interior  cell 
pressure,  to  the  ghost  cell,  Pg\ 


Pg  =  Pi  =P 


(3.76) 


To  set  the  temperature  of  the  wall  to  the  desired  temperature,  an  extrapolation 
from  inside  the  domain  to  outside  is  performed; 


T7 


(3.77) 


or 


7^  =  27^  —  7^ 

g  w  ^  i 


(3.78) 
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or,  in  terms  of  pressure  and  temperatures 


TP 

Pg 


2T 


TP 

Pi 


(3.79) 


Rearranging  for  we  obtain 

Pc  =  sfrr 

(7P) 

Immediately,  one  major  problem  can  be  seen  with  this  scheme,  in  that  a  value  of  wall 
temperature  of  1/2  the  adjacent  cell  temperature  will  cause  an  undehned  value,  or 
in  other  words,  this  expression  for  the  external  density  has  a  singularity.  Also,  lower 
temperatures  lead  to  an  unphysical  negative  density.  To  prevent  this  situation  in 
the  present  work,  the  wall  temperature  will  be  set  to  a  value  no  less  than  3/4  of  the 
temperature  in  the  adjacent  cell: 


T„. 


max 


(3.81) 


3.9  One  Dimensional  Shock  Solver 

The  first  test  of  bulk  viscosity  will  be  a  normal  shock  prohle,  using  a  one  dimen¬ 
sional,  direct  integration  of  the  Navier  Stokes  equations.  The  1-D  integration  code 
follows  the  method  outlined  in  Gilbarg  and  Paolucci  [37]. 

To  illustrate  the  effect  of  bulk  viscosity  on  shock  structure,  the  one-dimensional 
steady  state  Navier-Stokes  equations  are  solved  with  the  Runge-Kutta  ODE  solver. 
The  results  can  then  be  compared  to  the  observations  of  Alsmeyer  [25]  for  shock 
thickness  in  low-density  nitrogen.  The  procedure  of  Gilbarg  and  Paolucci  [37,59]  as 
implemented  by  Gamberos  and  Ghen  [60]  is  followed,  where  Navier-Stokes  equations 
in  one  dimension  are  solved: 


/ 


d 

dx 


pu 

pv?  +  p 


\ 


\  {Et+p)U  j 


d 

dx 


\p  -\-  pb) 


du 


\  {^y-l)Pr)  %  +  +  Pb)%U  j 


du  ^ 


(3.82) 
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which,  on  integrating,  become 


pu  =  A  (3.83) 

5  pT  /4  \  du  ^  /  .N 

^  f +  —  =  5  (3.84) 

/  ^  ,  a  dT  f  A  \  du  ^  ,  , 

(^‘  +  p>“-(7rT)75;*-(3^  +  f‘bf‘“*  =  ^  P.85) 

rearranging  these  for  =  du/dx  and  =  dT/dx,  we  obtain 

3  /  AT\ 

u.  =  --(b-A„-— )  (3,86) 

and 

~  +  ^AT  -C  +  Bu)  (3.87) 

p  \  2  4  y 

These  two  ODEs  are  integrated  using  an  adaptive  Runge-Kutta  scheme,  such  as  0DE45 
in  Matlab®or  NDSolve  in  Mathematica®.  The  latter  is  used  in  this  work. 

The  boundary  conditions  used  are  those  of  the  familiar  calorically  perfect  normal 
shock  relations  [41:86-94]: 


^  _  (7  -  1)M^ 

Pi  “  2  +  (7-1)M2) 

^  _  7-1  2 

Ml  7  -|-  1  (7  -|-  1)M2 


(3.88) 

(3.89) 


Note  that  continuity,  equation  (3.83),  requires  equation  (3.89)  to  be  the  reciprocal  of 
equation  (3.88). 


Pi  7-^1  7  -h  1 

_  h  _  ^fh 

Ti  hi  Pi  p2 


(3.90) 

(3.91) 


Note  that  relationships  for  a  more  general  case  than  the  calorically  perfect  gas  can  also 
be  written  in  terms  of  “thermodynamic  variables”  only,  as  shown  in  Anderson  [41:98- 
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101],  giving  the  commonly  known  Rankine-Hugoniot  equation: 


62  -  ei  = 


(3.92) 


where  v  is  the  specihc  volume.  This  equation  holds  no  matter  the  thermodynamic 
model,  however,  the  states  1  and  2  are  equilibrium  states. 

3.10  Axisymmetric  Test  Cases 

The  second  two  test  cases  involve  run  5  and  run  7  from  the  CUBRC  [1],  for 
which  experimental  wall  data  for  pressure  and  heating  is  available.  Solutions  are 
computed  with  AFIT-2D,  and  compared  to  experiment  and  DSMC  solutions  from 
Carr  [11].  The  effect  on  the  surface  pressure  and  heating  to  changes  in  bulk  viscosity 
will  be  analyzed.  Freestream  properties  are  listed  in  Table  2.  The  geometry  of  run  5 
is  depicted  in  Figure  6,  and  that  of  run  7  is  depicted  in  Figure  7. 


Table  2:  Freestream,  Wall  Conditions  for  CUBRC  Data 


Run  Ma  Re  (K)  Pqq  (Pa)  pop  (kg/m^)  Uqq  (m/s)  T^gu  (K)  Kn 

5  15.3  23746  52.28  2.523  0.00016  2252.47  303.33  0.0014 

7  15.6  26624  42.61  2.227  0.00018  2072.64  297.22  0.0014 


Figure  6:  Geometry  of  Run  5,  from  Holden  and  Wadhams  [1] 
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IV.  Results 


4..I  Shock  Structure 

The  inclusion  of  the  bulk  viscosity  model  of  Capelletti  et  al  in  the  Navier-Stokes 
equations  improves  shock  thickness,  as  can  be  seen  in  Figure  8.  This  hgure  shows 
inverse  shock  thickness,  as  dehned  by 


—  =  max 
t.. 


dpr, 


dx/ Al 


(4.1) 


where 


p„  =  (4.2) 

P2  -  Pi 

Here,  pi  is  the  pre-shock  density,  p2  is  the  post-shock  density,  and  Ai  is  the  pre-shock 
(or  freestream)  mean  free  path,  dehned  by  Alsmeyer  [25:499-500].  The  red  line  is 
the  inverse  shock  thickness  calculated  from  Alsmeyer’s  original  electron-beam  density 
observations.  These  observations  were  taken  at  p  =  50  mTorr,  or  6.67  Pa,  and  at 
a  temperature  of  300K.  Note  that  the  low  Mach  number  shocks  are  relatively  thick. 
The  shock  thickness  decreases  to  a  minimum  between  Ma  =  4  and  Ma  =  6,  and  then 
increases  again  as  the  Mach  number  increases.  The  DSMC  solution,  from  Bird,  agrees 
quite  well  with  these  measurements.  The  Navier-Stokes  solutions,  calculated  using  the 
Sutherland  model  for  shear  viscosity,  and  no  bulk  viscosity  predict  a  thinner  shock 
than  the  observations.  They  predict  a  decreasing  shock  thickness  with  increasing 
Mach  number,  through  the  Mach  number  range  of  the  observations.  A  large  error 
is  seen,  comparing  to  the  observations.  The  constant  bulk  viscosity  solutions,  while 
closer  to  the  observations,  still  predict  a  thinning  of  the  shock  throughout  the  observed 
range.  The  Navier-Stokes  solutions  with  the  temperature-dependent  bulk  viscosity 
model  from  Cappelletti  et  al  [36]  present  a  signihcant  improvement  over  the  previous 
two,  for  Mach  numbers  up  to  8.  At  Ma  =  9  and  Ma  =  10,  this  model  predicts  thinner 
shocks,  almost  the  same  thickness  as  those  of  the  constant  bulk  viscosity  solutions.  It 
should  be  noted  the  post-shock  temperatures  exceed  the  range  of  temperatures  of  the 
curveht  equation  2.24,  and  indeed,  far  exceed  the  temperature  at  which  vibrational 
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Figure  8:  Effect  of  bulk  viscosity  on  shock  thickness 

modes  are  expected  to  dominate  the  nonequilibrium  behavior  of  the  flow.  This  same 
problem  does  not  occur  with  the  two  dimensional  examples  discussed  below,  due  to 
their  very  low  freestream  temperature. 

Figures  9  (a)  and  (b)  show  the  normalized  density  profiles  as  defined  by  Eq  (4.2), 
at  Ma  =  4  and  Ma  =  6.  Here  again,  for  the  bulk  viscosity  models,  the  thicker 
shock  can  be  seen,  but  also  the  distance  of  the  shock  from  the  downstream  boundary 
condition  can  be  seen.  This  distance  can  be  though  of  as  a  “relaxation  distance.” 
The  relaxation  time  present  in  the  rotational  collision  number  is  physically  visible 
in  the  shock  solution.  The  two  are  related  by  convection.  It  is  worth  mentioning 
again  that  the  jump  conditions  have  not  changed,  including  those  for  entropy.  The 
thicker  shock  may  give  a  lower  peak  entropy  generation.  A  note  of  caution,  however, 
is  in  order.  Any  claim  of  the  validity  of  the  bulk  viscosity  model  based  on  continuum 
breakdown  parameters,  such  as  entropy  generation,  would  be  premature  speculation 
without  further  study. 
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(a)  Ma  =  4  (b)  Ma  =  6 

Figure  9:  Normalized  density  profiles 

4-2  Hollow  Cylinder  Flare  Analysis 

The  hollow  cylinder  flare  flow,  in  the  cylinder  region,  is  very  similar  to  a  flat 
plate  flowfield.  The  surface  of  the  body  is  aligned  with  the  flow.  The  researchers  at 
CUBRC  have  carefully  created  a  sharp  leading  edge,  sharp  enough  that  no  lip  radius 
is  given.  This  situation  results  in  a  nearly  attached  leading  edge  shock,  accompanied 
by  a  very  high  temperature  and  heating  rate.  This  portion  of  the  body  is  where  the 
constant  temperature  wall  listed  in  the  CUBRC  run  5  conditions,  shown  in  Table  2 
would  not  be  expected  to  hold.  In  fact,  none  of  the  thin-hlm  heat  transfer  gauges 
or  pressure  sensors  on  the  body  are  located  near  this  point  [1].  A  boundary  layer 
forms  on  the  surface,  with  a  weak  shock  bounding  its  upper  surface.  This  boundary 
layer  contains  high  temperatures,  as  the  flow  is  slowed  considerably  in  this  boundary 
layer  from  the  freestream,  and  the  kinetic  energy  contained  therein  is  converted  to 
thermal  and  internal  energy,  here  in  the  form  of  rotational  energy.  To  be  fair,  the 
region  near  the  stagnation  point  should  be  hot  enough  to  excite  some  vibration,  but 
its  overall  effect  on  the  flow  is  here  assumed  to  be  small.  The  density  predicted  in 
the  continuum  simulations  in  this  boundary  layer  is  actually  lower  than  that  of  the 
freestream.  This  feature  can  be  seen  in  Figure  11  (a)  through  (c).  In  the  calorically 
perfect  gas  assumption,  the  thermodynamic  pressure,  density,  and  temperature  are 
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related  by  T  oc  7p/p,  where  7  is  the  ratio  of  specihc  heats.  When  temperature 
increases,  at  nearly  the  same  pressure  (see  the  boundary  pressure  contour  plots  in 
Figure  10),  the  density  must  decrease. 
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(a)  fiB  =0 


(b)  fiB  =  0.8/i 


(c)  fj,B  from  Cappelletti  et  al  [36] 


(d)  DSMC  from  Carr  [11] 


Figure  10:  Hollow  cylinder  pressure 
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(c)  fjLB  from  Cappelletti  et  al  [36] 

Figure  11:  Hollow  cylinder  density 
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The  next  feature  present  in  the  continuum  solutions,  is  a  gradual  increase  in 
density  and  pressure  as  the  boundary  layer  is  turned  at  the  start  of  the  flare,  as 
can  be  seen  in  Figures  11  and  10.  The  leading  edge  of  this  pressure  and  density 
increase,  begins  to  exhibit  a  shock-like  structure,  when  following  it  toward  its  apparent 
intersection  with  the  boundary  layer  edge  shock.  This  feature  is  best  seen  in  the 
density  contour  plots  of  Figure  12.  This  region  is  an  example  of  shock-boundary  layer 
interaction.  The  problem  of  shock-boundary  layer  interaction  at  hypersonic  speed  has 
caused  serious  damage  to  vehicles.  For  example,  the  X-15,  while  testing  a  dummy 
ramjet  attached  a  pylon  on  the  bottom  of  the  aircraft,  experienced  a  loss  of  this  pylon 
due  the  shock  from  the  ramjet  model  impinging  on  the  pylon  boundary  layer  [17:562- 
565].  This  region  is  one  of  very  high  pressure  as  can  be  seen  in  Figure  10.  Figure  13 
depicts  pressure  coefficient  as  a  function  of  axial  distance  from  the  leading  edge  of  the 
flare.  The  CFD  predictions  of  surface  pressure  match  those  experimentally  observed 
at  CUBRC  quite  well  in  cylinder  region.  The  pressure  rise  on  the  flare,  through  the 
shock  impingement  point  is  slightly  underpredicted  by  all  three  CFD  solutions.  The 
portion  of  the  surface  aft  of  the  impingement  shows  an  underprediction  of  pressure  in 
the  expansion  region.  It  should  be  noted  that  this  feature  was  also  predicted  in  a  hner 
grid  with  AFIT-2D,  and  also  in  the  much  hner,  unstructed  grid  FLUENT ©solution  of 
Carr  [11:52].  From  these  observations,  it  can  be  speculated  that  the  disagreement  in 
pressure  in  this  expansion  region  is  less  related  to  grid  convergence,  and  more  related 
to  a  shortcoming  in  the  model  of  the  physics  of  the  CFD  and  DSMC  solutions. 
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(c)  fjLB  from  Cappelletti  et  al  [36] 

Figure  12:  Hollow  cylinder  density  contours 
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Figure  13:  Hollow  cylinder  surface  pressure  coefficient 

Figure  14  depicts  Stanton  number,  or  normalized  heat  transfer  at  the  wall. 
In  this  hgure  can  be  seen  that  the  CFD  predictions  of  heat  transfer  are  not  far 
the  experimentally  observed.  The  bulk  viscosity  models  both  seem  to  provide  some 
improvement  in  heat  transfer  prediction  over  the  zero  bulk  viscosity  CFD,  in  the 
forward  side  of  the  imingement  point,  and  in  heat  transfer  predicted  in  aft  of  the 
impingment  point,  near  the  end  of  the  body. 
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Figure  14:  Hollow  cylinder  Stanton  number 
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Figure  15:  Hollow  cylinder  velocity  vector  detail:  /i^  =  0 

The  downstream  side  of  the  interaction  point  includes  two  features.  The  hrst 
being  the  outer  edge  of  the  shock  layer,  which  appears  to  approach  the  slope  of  the 
shock  developing  out  of  the  boundary  layer,  mentioned  above.  The  second  feature  is 
less  like  a  shock.  The  momentum  and  energy  flowing  through  this  area  form  a  conical 
“sheet”  of  high-velocity  flow,  in  the  expanding  region  between  the  flare  and  the  edge 
of  the  shock  layer.  The  velocity  vectors  in  the  closeup  of  hgure  15  show  the  detail 
of  this  feature.  Traversing  the  flare  portion  of  the  shock  layer  edge,  a  higher  than 
freestream  density  part  of  the  flow  can  be  seen  at  the  interaction  point.  This  region 
is  the  cylinder  portion  of  the  shock  layer.  A  gradual  relaxing  of  the  shock  angle  can 
be  seen  moving  downstream.  This  high  pressure  high  density  region  expands  and 
accelerates  between  the  shock  layer  edge  and  the  wall.  The  effect  of  bulk  viscosity  on 
these  solutions  is  visible  in  Figures  12  and  16.  Figures  12  (b)  and  (c)  depict  a  wider 
shock  at  the  edge  of  the  shock  layer  than  does  Figure  12  (a).  Also,  the  shock  feature 
developing  in  the  boundary  layer,  ahead  of  the  interaction  point  is  thicker.  Overall 
shock  layer  thickness  is  not  greatly  affected. 
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(a)  AiB  =  0  (b)  ^JLB  =  0.8/i 


(c)  fjLB  from  Cappelletti  et  al  [36] 

Figure  16:  Hollow  cylinder  pressure  contours 
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Figure  17:  Hollow  cylinder  pressure  comparison  locations 

Figure  17  depicts  several  grid  line  cuts,  where  pressure  information  has  been 
extracted  from  both  the  CFD  and  the  DSMC.  The  hrst  cut  is  taken  ahead  of  the 
compression  corner.  The  second  is  taken  just  after  the  corner,  where  there  should  be 
higher  pressure  at  the  wall.  The  third  cut  is  taken  where  the  compression  should  be 
well-developed,  just  before  the  high  pressure  shock- impingement  region.  The  fourth 
cut  is  near  the  maximum  pressure  point.  The  hfth  cut  is  taken  past  the  impingement, 
capturing  the  high-velocity  “sheet.” 
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Figure  18:  Extraction  of  pressure  data,  run  5,  x/L  =  0.32,  i  =  175 

Figures  18  through  22  represent  extractions  of  pressure  contours  taken  on  con¬ 
stant  values  of  or  i,  shown  in  Figure  17,  from  the  three  CFD  solutions,  and  the 
DSMC  solutions  from  Carr  [11],  Carr’s  DSMC  solutions  here  presented  were  com¬ 
puted  with  Boyd’s  MONACO  DSMC  solver  [61].  Figure  18  compares  pressure  prohles 
between  the  three  CFD  solutions,  and  the  MONACO  DSMC  solution  computed  by 
Carr  [11].  This  extraction  of  pressure  information  was  taken  in  the  cylinder  region, 
slightly  ahead  of  the  flare.  The  CFD  solutions  predict  a  slightly  thicker  shock  layer, 
and  an  overall  thinner  shock  than  does  MONACO.  The  CFD  solutions  computed  with 
both  constant  bulk  viscosity  and  the  temperature  dependent  bulk  viscosity  predict 
a  slightly  thicker  shock,  and  a  slightly  thinner  shock  layer  than  does  the  zero  bulk 
viscosity  CFD  solution.  The  MONACO  solution  predicts  a  local  pressure  peak  near 
the  wall,  a  feature  not  predicted  by  the  CFD  solutions.  Figure  19  compares  pressure 
prohles  along  constant  grid  coordinate  i  =  263.  This  extraction  of  pressure  informa- 
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Figure  19:  Extraction  of  pressure  data,  run  5,  x/L  =  0.48,  i  =  263 

tion  was  taken  near  the  junction  of  the  cylinder  and  flare.  The  CFD  solutions  predict 
a  slightly  thicker  shock  layer,  and  an  overall  thinner  shock  than  does  MONACO.  The 
shock  thickness  in  the  case  of  both  constant  bulk  viscosity  and  that  of  the  temperature 
dependent  model  predict  a  slightly  thicker  shock,  and  a  slightly  thinner  shock  layer 
than  .  The  MONACO  solution  predicts  slightly  higher  pressure  at  the  wall  than  that 
predicted  by  the  three  CFD  solutions.  There  is  a  qualitative  similarity  of  the  four 
prohles  not  seen  in  the  i  =  175  cut  of  Figure  18.  Figure  20  compares  pressure  prohles 
at  i  =  348,  or  slightly  ahead  of  the  shock-boundary  layer  interaction  region,  as  seen  in 
Figure  17.  The  three  CFD  solutions  predict  a  thinner  shock  than  does  the  MONACO 
DSMC,  and  here  a  slightly  thicker  shock  layer.  The  presence  of  two  shocks  can  be 
seen  here.  The  outer,  weak  boundary  layer  shock  of  the  cylinder  region  can  be  seen  in 
the  upper  left.  Note  that  the  constant  bulk  viscosity  and  temperature  dependent  bulk 
viscosity  have  brought  this  shock,  and  hence  the  overall  shock  layer  thickness  down 
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Figure  20:  Extraction  of  pressure  data,  run  5,  x/L  =  0.62,  i  =  348 

closer  to  that  of  the  MONACO  DSMC.  This  shock  thickness  is  also  slightly  increased 
by  the  inclusion  of  bulk  viscosity.  The  inner,  stronger  shock  is  predicted  closer  to 
the  body  than  the  DSMC  solution  by  all  three  CFD  solutions.  The  inner  shock  is 
also  predicted  to  be  thinner  than  the  inner  shock  predicted  by  DSMC.  The  effect  of 
bulk  viscosity  on  this  inner  shocks  thickness  is  minimal.  The  wall  pressure  predicted 
by  CFD  solutions  including  bulk  viscosity  is  higher  than  without,  and  closer  to  that 
predicted  by  DSMC.  Figure  21  compares  pressure  prohles  in  the  shock  impingement 
region,  where  pressure  is  expected  to  be  highest,  as  depicted  in  Figure  17.  The  three 
CFD  pressure  prohles  again  show  a  similar  shape  to  that  of  the  DSMC.  The  shock 
predicted  by  DSMC  is  thicker,  as  is  the  shock  layer.  The  constant  and  temperature- 
dependent  bulk  viscosity  CFD  solutions  predict  a  thicker  shock  than  does  the  CFD 
solution  computed  with  no  bulk  viscosity.  The  pressure  at  the  surface  is  predicted  by 
all  three  CFD  solutions  is  higher  than  that  predicted  by  DSMC.  It  should  be  noted 


Figure  21:  Extraction  of  pressure  data,  run  5,  x/L  =  0.67,  i  =  380 

that  the  CFD  surface  pressure  is  in  closer  agreement  to  the  experimentally  observed 
surface  pressure  [11:52],  Figure  22  compares  pressure  prohles  past  the  impingement 
point.  This  extraction  is  represented  as  the  aft-most  white  line  in  Figure  17.  It  should 
be  mentioned  that  at  this  point  the  CFD  solution  is  in  better  agreement  with  exper¬ 
imental  surface  pressure  than  is  the  DSMC  solution  of  Carr,  and  the  DSMC  solution 
predicts  the  shock-boundary  layer  interaction  point  to  be  farther  upstream  than  do 
the  CFD  solutions  [11:52].  The  CFD  solutions  predict  a  slightly  thinner  shock  and 
shock  layer  than  does  the  DMSC.  The  effect  of  bulk  viscosity  in  the  CFD  solutions  is 
a  thicker  predicted  shock  than  that  predicted  by  the  zero  bulk  viscosity  CFD  solution. 
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Figure  22:  Extraction  of  pressure  data,  run  5,  x/L  =  0.83,  i  =  450 
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4-5  Double  Cone  Analysis 

The  double  cone  flow  is  slightly  more  complicated  than  that  of  the  hollow  cylin¬ 
der  flare.  A  thin  shock  layer  is  present  on  the  forward  25°  cone.  As  in  the  previous 
solution,  the  leading  edge  is  sharp.  It  should  be  noted  again  that  the  continuum  so¬ 
lutions  presented  here  have  been  calculated  on  a  model  modified  to  have  a  very  small 
nose  radius.  This  choice  was  made  to  avoid  creating  a  grid- singularity  boundary  con¬ 
dition,  and  modeling  the  nose  as  sharp.  The  nose  radius  used  was  0.0625  inches. 
Approaching  the  corner  to  the  55°  rear  cone,  an  increase  in  pressure  and  density  is 
seen  ahead  of  the  corner. 

The  shock  that  develops  in  this  region  interacts  with  the  boundary  layer  flow, 
separating  the  flow  into  a  laminar  separation  region.  Figure  23  depicts  this  region 
with  a  velocity  vector  presentation.  The  blue  lines  in  this  figure  represent  the  bounds 
of  snbsonic  flow  within  the  solution.  As  with  most  shock/bonndary  layer  interactions, 
a  clearly  defined  shock  is  not  present  near  the  wall,  as  can  be  seen  in  Figure  25.  This 
shock  feature  is  considerably  thicker  than  that  seen  upstream. 

Between  the  shock  layer  edge  and  the  separated  region,  the  mass  and  momen¬ 
tum  of  the  forward  shock  layer  flows.  The  forecone  shock  smoothly  merges  into  the 
“separation  shock,”  which  bonnds  this  region.  The  separated  region  can  be  seen  in 
Figure  23 

Near  the  end  of  the  separated  region,  the  shock  layer  edge  splits.  One  branch 
becomes  a  nearly  normal  shock,  ahead  of  the  rear  cone,  and  the  other  connects  to 
the  reattachment  of  the  separated  region.  Nompelis  et  al  [62]  call  this  the  “trans¬ 
mitted”  shock.  The  point  of  maximnm  pressure  is  just  aft  of  this  transmitted  shock. 
Downstream  of  this  point,  the  maximnm  pressure  region  produces  a  “sheet”  of  high 
velocity,  described  by  Nompelis  et  al  as  a  jet.  This  sheet  of  high  speed  flow  is  much 
closer  to  the  wall  than  that  fonnd  in  the  hollow  cylinder  flare  interaction.  The  upper 
surface  of  this  jet  or  sheet  is  a  contact  snrface. 
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Figure  23:  Double  cone  separation  region  detail 
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(a)  /is  =  0 


(b)  ^jlb  =  0.8/i 


(c)  fj,B  from  Cappelletti  et  al  [36] 

Figure  24:  Double  cone  density 

Figure  24  compares  contours  of  density.  The  shock  layer  is  progressively  thicker 
on  moving  from  the  solution  with  no  bulk  viscosity,  to  the  that  of  constant  bulk 
viscosity  and  is  thickest  in  the  temperature-dependent  model. 
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(a)  /is  =  0 


(b)  ^jlb  =  0.8/i 


(c)  fj,B  from  Cappelletti  et  al  [36] 

Figure  25:  Double  cone  density  contours 

The  shock  thickness  is  most  visible  in  Figure  25.  The  shock,  and  shock  layer 
increases  in  thickness  when  moving  from  the  zero  velocity  solution  to  that  of  the  con¬ 
stant  bulk  viscosity,  and  hnally  to  the  temperature  dependent  bulk  viscosity  solution. 
The  separation  region  forward  of  the  corner  can  be  seen  in  the  density  contour  plot. 
The  constant  bulk  viscosity  solution  predicts  a  larger  separation  region,  extending 
slightly  forward  of  that  predicted  by  the  zero  bulk  viscosity  solution.  The  temper¬ 
ature  dependent  bulk  viscosity  solution  extends  even  farther  forward,  and  is  larger. 
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(a)  /is  =  0 


(b)  ^jlb  =  0.8/i 


(c)  fj,B  from  Cappelletti  et  al  [36] 
Figure  26:  Double  cone  pressure 


Figure  26  depicts  the  pressure  field  of  the  three  CFD  solutions.  The  separated 
flow  region  can  be  seen,  as  a  region  of  increased  pressure  compared  to  the  shock  layer 
on  the  25°  forecone.  The  separated  region  is  predicted  to  be  larger  in  the  constant 
bulk  viscosity  solution  than  in  the  zero  bulk  viscosity  solution.  The  largest  separation 
region  can  be  seen  in  the  temperature  dependent  bulk  viscosity  solution.  The  latter 
solution  also  predicts  a  higher  maximum  pressure,  slightly  aft  of  maximum  pressure 
predicted  by  the  zero  and  constant  bulk  viscosity  solutions. 
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(c)  fj,B  from  Cappelletti  et  al  [36] 

Figure  27:  Double  cone  pressure  contours 


Figure  27  depicts  pressure  contours  for  the  three  CFD  solutions.  The  shock 
layer  thickness  is  greater  in  the  temperature  dependent  bulk  viscosity  solution,  than 
in  the  other  two  CFD  solutions.  The  larger  separated  flow  region  near  the  corner 
can  be  seen  in  the  nonzero  bulk  viscosity  solutions,  the  largest  seen  in  the  case  of  the 
temperature  dependent  bulk  viscosity.  The  highest  peak  pressure  can  also  be  seen 
clearly  in  Figure  27  (c). 
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Figure  28:  Double  cone  surface  pressure  coefficient 


The  pressure  coefficient  on  the  surface  is  depicted  in  Figure  28.  Note,  the 
separation  region  is  slightly  better  predicted  by  the  constant  bulk  viscosity  model. 
The  solution  computed  with  temperature  dependent  bulk  viscosity  predicts  slightly 
earlier  separation,  but  still  separates  aft  of  the  experimental  measurements,  as  do  the 
other  CFD  solutions,  and  the  DSMC  solutions  of  Carr  [11:64].  The  three  pressures 
near  and  aft  oi  x/L  =  0.4  are  predicted  by  the  temperature  dependent  bulk  viscosity 
in  near  perfect  agreement  with  the  observed  pressures,  as  are  the  two  pressures  on 
the  forward,  rising  pressure  side  of  the  impingement  point  near  x/L  =  0.5.  Aft  of 
the  impingment  point,  the  three  CFD  solutions  exhibit  more  oscillation  in  pressure 
than  that  seen  in  the  DSMC  solutions.  There  appears  to  be  some  scatter  in  the 
experimental  results  in  this  region.  The  Stanton  number  is  depicted  in  Figure  29. 
The  three  CFD  solutions  are  depicted  with  dashed  lines,  to  show  that  the  associated 
predictions  of  heating  are  nearly  identical.  The  heat  transfer  is  underpredicted  on 
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Figure  29:  Double  cone  Stanton  number 

the  forecone.  The  heat  transfer  predicted  by  CFD  through  the  separation  region, 
from  x/L  =  0.4  to  x/L  =  0.5  is  in  agreement  with  the  observed  heat  transfer.  The 
two  heating  predicitons  just  past  x/L  =  0.5  are  in  agreemement  with  the  observed 
heat  transfer.  Unfortunately,  no  experimental  point  is  located  closer  to  the  predicted 
maximum  heat  transfer.  The  heat  transfer  on  the  portion  of  the  55°  cone  aft  of  the 
impingement  point  is  underpredicted  by  the  three  CFD  solutions.  Figure  30  depicts 
locations  of  pressure  cuts  taken  to  compare  the  CFD  solutions  to  the  DSMC  solutions 
of  Carr  [11:64-65].  For  the  double  cone  case,  solutions  from  both  MONACO,  which 
uses  the  Parker  rotational  relaxation  model,  and  DS2V,  Bird’s  solver  (of  [9]),  using  a 
constant  rotational  collision  number,  are  available.  Carr  reports  the  latter  solutions 
to  be  in  closer  agreement  to  experiment  than  the  former  [47].  Figures  31  through  36 
contain  selected  pressure  extractions  at  the  locations  depicted  in  Figure  30. 
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Figure  30:  Double  cone  pressure  comparison  locations 

Figure  31  depicts  extractions  of  pressure  data  from  the  three  CFD  solutions, 
MONACO  and  DS2V  DSMC  codes.  This  extraction,  as  seen  in  Figure  30  is  taken 
on  the  25°  forecone,  ahead  of  the  separation  region.  Both  DSMC  solvers  predict  a 
thicker  shock  layer,  and  thicker  shock  than  that  seen  in  the  CFD  solutions.  The  CFD 
solutions  computed  with  constant  and  temperature  dependent  bulk  viscosity  predict  a 
thicker  shock  than  does  the  CFD  solution  computed  with  no  bulk  viscosity.  Increased 
pressure  at  the  wall  is  visible  in  the  DS2V  solution  only.  Note  that  the  two  DSMC 
solutions  predict  different  pressure  at  this  point. 

Figure  31  depicts  extractions  of  pressure  data  taken  from  near  the  forward  end 
of  the  separation  region.  At  this  location,  x/L  =  0.4,  the  CFD  solution  computed- 
with  temperature  dependent  bulk  viscosity  is  the  only  solution  which  agrees  with  the 
observed  pressure,  as  can  be  seen  in  Figure  28.  The  CFD  solutions  with  bulk  viscosity 
predict  a  thicker  shock  than  does  the  CFD  solution  with  no  bulk  viscosity.  The  pres¬ 
sure  rise  toward  the  surface  is  seen  in  the  constant  bulk  viscosity  CFD  solution  and 
the  DS2V  DSMC  solution,  to  a  lesser  extent.  Recall  that  the  constant  bulk  viscosity 
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Figure  31:  Extraction  of  pressure  data,  run  7,  x/L  =  0.31,  i  =  169 
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Figure  32:  Extraction  of  pressure  data,  run  7,  x/L  =  0.40,  i  =  215 

and  a  constant  rotational  collision  number  are  similar  phenomenological  models  for 
rotational  relaxation,  so  some  similarity  is  expected.  The  hot  side  of  the  shock  is 
predicted  to  be  thicker  in  DS2V  than  with  the  CFD  solutions. 

Figure  33  depicts  pressure  extractions  taken  at  x/L  =  0.46,  within  the  separa¬ 
tion  region.  As  can  be  seen  in  Figure  28,  the  CFD  solution  computed  with  constant 
bulk  viscosity  nearly  agrees  with  the  observed  surface  pressure.  The  constant  and 
temperature  dependent  bulk  viscosity  model  and  DS2V  predict  one  shock,  whereas 
the  zero  bulk  viscosity  CFD  solution  and  MONACO  solutions  predict  the  forecone 
shockand  separation  shock  as  separate  entities  at  this  cut  location. 

Figure  34  depicts  pressure  extractions  from  x/L  =  0.52,  or  the  forward  side 
of  the  shock  impingement  point.  The  CFD  solution  computed  with  temperature 
dependent  bulk  viscosity  predicts  the  thickest  shock  layer,  and  thicker  shocks  than 
the  other  two  CFD  solutions.  The  thinnest  shock  layers  are  predicted  by  the  two 
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Figure  33:  Extraction  of  pressure  data,  run  7,  x/L  =  0.46,  i  =  250 
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Figure  34:  Extraction  of  pressure  data,  run  7,  x/L  =  0.52,  i  =  292 

DSMC  solutions.  Note  that  DS2V  predicts  a  pressure  maximum  above  the  surface  at 
this  location.  A  much  smaller  local  pressure  maximum  at  this  point  can  be  seen  in  the 
zero  and  constant  bulk  viscosity  CFD  solutions.  There  is  some  similarity  of  pressure 
prohles  seen  when  comparing  temperature  depdendent  bulk  viscosity  and  MONACO, 
which  uses  a  temperature  dependent  rotational  collision  number. 

Figure  35  depicts  pressure  prohles  slightly  downstream,  where  the  temperature 
dependent  model  overpredicts  presssure,  when  compared  to  the  observations  of  Holden 
and  Wadhams  [1].  Here  this  solution  again  predicts  a  thicker  shock  layer  than  do  the 
other  two  CFD  solutions.  The  zero  bulk  viscosity  CFD  solution  nearly  matches  the 
shock  location  predicted  by  the  MONACO  solution,  while  DS2V  predics  a  thinner 
shock  layer.  The  DS2V  and  zero  bulk  viscosity  CFD  solutions  predict  a  local  pressure 
maximum  at  this  cut  location. 
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Figure  35:  Extraction  of  pressure  data,  run  7,  x/L  =  0.53,  i  =  299 
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Figure  36:  Extraction  of  pressure  data,  run  7,  x/L  =  0.58,  i  =  344 


Figure  36  depicts  the  pressure  profiles  ai  x/L  =  0.58,  past  the  shock  impinge¬ 
ment  point.  The  DSMC  solutions  predict  the  thinnest  shock  layer  .  The  thickest  shock 
and  shock  layer  is  predicted  by  the  temperature  dependent  bulk  viscosity  model.  A 
thinner  shock  layer  is  predicted  by  the  constant  bulk  viscosity  model,  and  the  thinnest 
of  the  three  CFD  solutions  was  computed  with  zero  bulk  viscosity.  The  zero  and 
constant  bulk  viscosity  CFD  solutions  exhibit  some  similarities  to  the  DS2V  DSMC 
solution. 
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V.  Conclusions 


In  studying  high  speed  and  high  altitude  flows,  a  means  for  accounting  for  non¬ 
equilibrium  effects,  such  as  the  equilibration  process  between  internal  and  transla¬ 
tional  modes  of  energy  is  desired.  The  geometry  of  the  flow  features  in  non-equilibrium 
situations  is  affected  by  this  equilibration  process.  The  shock  layer  thickness  in  high 
speed  flight  is  affected  by  this  process.  Location  of  shock/boundary  layer  interaction 
points  is  also  affected  by  this  process. 

The  bulk  viscosity,  incorporated  into  CFD  models,  to  some  extent,  accounts  for 
the  effect  of  the  finite  time  of  the  rotation-translation  energy  transfer  on  the  flow. 
There  is  a  limit  to  how  far  from  equilibrium  bulk  viscosity  can  be  used  to  model  this 
process. 

The  bulk  viscosity  has  been  shown  here  to  improve,  in  some  circumstances,  the 
shock  thickness  and  prediction  of  heat  transfer  on  the  body.  Of  the  two-dimensional 
cases  studied,  the  surface  heat  transfer  distribution  of  the  hollow  cylinder  flare  shows 
improvement  on  including  the  bulk  viscosity  model  from  Cappelletti  et  al.  The  double 
cone  case  shows  some  improvement  in  the  prediction  of  the  onset  of  flow  separation 
at  the  corner  between  front  and  rear  cones. 

The  peak  pressure  tends  to  be  slightly  over-predicted  when  including  bulk  vis¬ 
cosity.  The  behavior  in  the  post-shock-impingement  region  of  the  double  cone  case 
appears  to  be  oscillatory,  and  this  behavior  is  not  in  agreement  with  the  experi¬ 
mental  data.  The  pressure  on  the  cylinder  flare  case  aft  of  the  impingement  region 
“undershoots”  the  measured  values,  and  thereafter  remains  lower.  These  two  ex¬ 
amples  suggest  that  the  bulk  viscosity  model  is  more  suited  to  compressions,  and 
behaves  differently  in  expansion  than  compression.  This  result  is  not  too  surpris¬ 
ing,  given  the  disagreements  in  the  literature  over  which  rotational  relaxation  models 
are  appropriate  for  a  particular  situation.  The  temperature  dependent  models  here 
presented  use  “sudden”  approximations,  which  assume  a  collision  duration  shorter 
than  the  rotational  period  of  the  molecule.  Bran  and  Jonkman  [30]  and  Wysong  and 
Wadsworth  [49]  point  out  that  if  the  opposite  is  assumed,  the  gas  does  not  behave 
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as  with  a  “single  relaxation  time.”  McCourt  et  al  [2:496]  describe  the  difficulties 
in  measuring  rotational  relaxation  with  ultrasonic  means,  and  the  difficulty,  at  least 
in  1990,  of  determining  the  rotational  temperature  using  light  scattering  methods, 
such  as  depolarized  Rayleigh  scattering.  Successful  determination  of  the  behavior  of 
rotational  energy  transfer  of  H2  has  been  performed  with  light  scattering,  but  it  is 
reported  that  the  greater  number  of  rotational  states  and  their  tighter  distribution 
for  heavier  molecules  such  as  N2  has  made  measurements  more  difficult.  It  is  hoped 
that  these  methods  will  continue  to  improve,  and  be  able  to  provide  more  information 
about  rotational  relaxation. 

It  has  been  remarked  that  the  bulk  viscosity  is  appropriate  for  “small”  deviations 
from  equilibrium  [63].  In  the  axisymmetric  solutions  presented  in  this  work,  that 
assessment  is  supported.  The  weaker-shock  run  5,  the  hollow  cylinder  flare,  shows  an 
improvement  in  agreement  with  observations  when  calculated  with  the  bulk  viscosity. 
The  stronger-shock  run  7,  the  double  cone,  shows  less  improvement.  The  Jeans  or 
Jeans-Landau- Teller  model  is  reported  to  be  more  appropriate  for  larger  deviations 
from  equilibrium  [44]. 

The  inclusion  of  the  internal  energy  modes  relaxation  in  determining  the  heat 
conduction  would  be  a  natural  next  step  for  future  work.  The  shear  viscosity  and 
heat  transfer  cross  sections  are  also  calculated  in  the  work  of  Cappelletti  et  al  [36]. 
The  effect  of  including  this  enhanced  heat  transfer  treatment,  and  shear  viscosity  is 
expected  to  be  small  in  the  case  of  the  conditions  and  geometries  studied  in  this  work. 
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In  the  work  of  Heck  and  Dickinson  [24],  on  which  the  work  of  Cappelletti  et 
al  [36],  it  is  mentioned  that  the  transport  coefficients,  such  as  shear  viscosity,  bulk 
viscosity,  and  heat  conduction  were  originally  meant  for  a  second  order  Chapman- 
Enskog  expansion  of  a  generalized  Boltzmann  Equation,  known  as  the  Curtiss-Kagan- 
Maksimov  equation.  Future  work  in  this  area  should  discover  the  form  of  these  second 
order  hydrodynamic  equations,  and  determine  how  these  equations  behave,  and  if  they 
share  the  same  problems  reported  for  their  second-order  hydrodynamics  counterpart, 
the  Burnett  equations,  obtained  from  a  similar  expansion  process  on  the  Boltzmann 
equation. 
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Appendix  A.  Align  Shock  Modification 


!  Added  7  Oct  07  Capt  Abram  Claycomb 

I  AT  1\T  O XJ n 7~* X<^  I 

I  j^  ^  ^yj  1^  yj  y^  i 

!  Aligns  Grid  to  and  clusters  points  around  shock,  keeping  adequate  boundary 
!  layer  resolution.  Taken  from  Gnoffo,  Hartung  and  Greendyke(AlAA-93-270-864) 


subroutine  align_shock() 

use  runtime_module 
implicit  none 


integer 

:  :  i,  j  ,  jt 

integer 

: :  j_star,  j_max 

integer 

: :  ncells_i,  nfaces_i 

integer 

: :  ncells_j ,  nfaces_j 

integer 

:  :  cell_imin,  cell_imax 

integer 

:  :  cell_jmin,  cell_jmax 

!  Index  counters 

!  j-location  of  end  of  b.l.  cells 


!  node  locations 

real (kind=RKlND) ,  dimension( : , : ) ,  pointer  : :  x,  y 
!  new  face  locations 

real (kind=RKlND) ,  dimension( : , : ) ,  allocatable  ::  xfnew,  yfnew 
!  cell  centers 

real (kind=RKlND) ,  dimension( : , : ) ,  allocatable  ::  xf_c,  yf_c  !,  x_c,  y_c 
!  change  in  x,  y  from  face  to  face  in  j  direction 
real (kind=RKlND) ,  dimension( : , : ) ,  allocatable  ::  dnx,  dny 
!  change  in,  'curvilinear'  normal  distance  in  j  away  from  wall 
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real (kind=RKIND) ,  dimension( : , : ) ,  allocatable  ::  dnf,  nold,  noldf 

real (kind=RKIND) ,  dimension( : , : ) ,  allocatable  ::  dnhat,  nhat,  anew,  nnewf 

!  running  total  of  distance  to  last  cell  center,  outer  face  of  grid 
real (kind=RKlND) ,  dimension( : ) ,  allocatable  ::  nlastcell,  nlastface 

!  shock  location  array 

real (kind=RKlND) ,  dimens i on( :) ,  allocatable  ::  nshk 


!  constants  for  the  code 
real (kind=RKlND)  ::  pi,  dpi 


!  shock,  b.l.  fitting  parameters 

real (kind=RKlND)  ::  Cl,  C2,  C2_max,  hminl,  hmin2,  rhowall,  asnd,  ep 
real (kind=RKlND)  : :  nxyf actor 

!  grid  metric  computation. . .  and  j  max  node  computation. . . 
real (kind=RKlND)  ::  dxl,  dyl,  dx2,  dy2,  na 


pi  =  acos(-ONE) 

ncells_i  =  num_i_cells (grid)  ! 
nfaces_i  =  num_i_f aces (grid)  ! 
cell.imin  =  CELL_MIN_1NDEX  ! 
cell_imax  =  cell_max_i_index(grid)  ! 

ncells_j  =  num_j_cells(grid) 
nfaces_j  =  num_j_f aces (grid) 
cell.jmin  =  CELL_MIN_INDEX 


internal  cells  only 

faces  =  #  cells  +  n_ghost 

includes  ghost  cells  (0  or  1  -  nghost) 

includes  ghost  cells  (ncells_i  +  nghost) 
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cell_jmax  =  cell_max_j_index(grid) 

allocateC  xf _c(l :ncells_i , 1 :nf aces_j )  ) 
allocateC  yf _c(l :ncells_i , 1 :nf aces_j )  ) 
allocateC  xfnewCl :ncells_i , 1 :nf aces_j )  ) 
allocateC  yfnewCl :ncells_i , 1 :nf aces_j )  ) 

!  allocateC  x_cCl :ncells_i , 1 :ncells_j )  ) 

!  allocateC  y_cCl :ncells_i , 1 :ncells_j )  ) 
allocateC  dnxCl :ncells_i, 1 :ncells_j)  ) 
allocateC  dnyCl :ncells_i, 1 :ncells_j)  ) 
allocateC  dnf Cl :ncells_i, 1 :ncells_j)  ) 

allocateC  noldCl :ncells_i , 1 :ncells_j )  )  !  diet  from  wall  to  this  cell  center 

allocateC  noldf Cl :ncells_i , 1 :nf aces_j )  )  !  ...  face 

allocateC  dnhat Cl :ncells_i , 1 :ncells_j )  ) 

allocateC  nhat Cl :ncells_i , 1 :nf aces_j )  ) 

allocateC  nnewCl :ncells_i , 1 :ncells_j )  ) 

allocateC  nnewf Cl :ncells_i , 1 :nf aces_j )  ) 

allocateC  nlastcellCl :ncells_i)  ) 

allocateC  nlastf aceCl :ncells_i)  ) 

allocateC  nshkCl :ncells_i)  ) 

X  =>  grid'/oX 
y  =>  grid'/.y 

!  get  face  center  position  of  j-normal  faces... 

xf_c  =  HALF*C  xC2 :nf aces_i , 1 :nf aces_j )  +  xCl :ncells_i , 1 :nfaces_j )  ) 

yf_c  =  HALF*C  y C2 : nf aces_i , 1 : nf aces_j )  +  yCl :ncells_i , 1 :nfaces_j )  ) 

!  get  cell  centers... 
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!  x_c  =  HALF*(  xf_c(l :ncells_i,2:nfaces_j)  +  xf _c(l :ncells_i , 1 :ncells_j )  ) 
!  y_c  =  HALF*(  yf_c(l :ncells_i,2:nfaces_j)  +  yf _c(l :ncells_i , 1 :ncells_j )  ) 

!  get  distance  between  j-face  centers... 

dnx  =  xf _c(l :ncells_i , 2 :nf aces_j )  -  xf _c (1 :ncells_i , 1 :ncells_j ) 
dny  =  yf _c(l :ncells_i , 2 :nf aces_j)  -  yf _c (1 :ncells_i , 1 :ncells_j ) 

dnf  =  sqrt(dnx**2  +  dny**2) 

nlastcelK:)  =  HALF*dnf  ( :  ,  1) 
nlastface(:)  =  dnf(:,l) 
noldf (: ,1)  =  ZERO 

nold(:,l)  =  nlastcelK:) 
noldf(:,2)  =  nlastface(:) 

!  get  cumulative  old  curvilinear  distance... 
do  j  =2,  ncells_j 

!jt  =  min(j-l,l) 

!nold(:,j)  =  HALF*  (dnf  (:  ,jt)+dnf  (:  ,j ) )  +  nlastcelK:) 
nold(:,j)  =  HALF*dnf ( : , j)  +  nlastface(:) 
nlastcelK:)  =  nold(:,j) 

noldf (:,j+l)  =  dnf(:,j)  +  nlastface(:) 
nlastface(:)  =  noldf (:,j+l) 

end  do 

!!$  open(unit=101 ,  file='ntest .dat^ ,  status  =  'replace’) 
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!!$  do  j  =  1,  ncells_j 

!!$  writedOl,  *)  'noldf  (1 ,  d  j"'"!  d  )  =  ’>  noldf(l,  j+D 

!!$  writeClOl,  *)  'nold(l,  ’  >  j  > ’)  =  nold(l,j).  ’  ,rho  =  gridy„Q(R_,  1 » j ) 

!  !  $  end  do 

!!$  ! write (*,  ’(A)’)  'outer  face  distance  test  complete' 

!  density  hard-coded,  could  easily  add  runtime_module  variable  for  this, 

!  initial_module  input  line  for  input  file  of  which  var.  to  sense... 
call  shock_loc (nold,  nshk,  R_)  !  locate  s  value  of  shock  throughout  grid. . . 

Ido  i  =  1,  ncells_i 

!  write(*,  '(A,  I,  A,  F8.3)')  'i  ',  i,  '  nshk  =  ',  nshk(i) 

! end  do 

! write (*,  '(A)')  'Shock  Location  Test  Completed' 

j_star  =  f_star  *  ncells_j 

dpi  =  pi  /  j_star 

nhat(:,l)  =  ZERO  !  nondimensional  distance  of  wall  face  from  wall 
nnewf(:,l)  =  ZERO  !  dimensional  distance  of  wall  face  from  wall  (itself...) 

do  i  =  ncells_i,  1,  -1 


rhowall  =  grid’/oOCR.,  i ,  1) 


asnd  =  sqrt  (GAMMA*gr  idZQ  (P_ ,  i ,  1 )  /gr  idy„Q  (R_ ,  i ,  1 ) ) 


hminl  =  Re_cellw  *  grid'/oinuCi ,  1)  /  (rhowall  *  asnd  *  noldf  (i ,nf aces_j )  ) 
! write (*,*)  i,  Re_cellw,  gridymuCi , 1) ,  hminl 

hmin2  =  ONE  /  nfaces_j 

hminl  =  min (hminl,  hmin2) 

dnhat(i,l)  =  hminl 

nhat(i,2)  =  dnhat(i,l) 

Cl  =  (  f_star  /  dnhat(i,l)  )  **  (  ONE  /  j_star  )  -  ONE 

!write(*,*)  i,  Cl 

C2  =  max(  ONE,  (ONE  +  Cl*sin(dpi))  ) 

C2_max  =  C2 

!dnhat(i,2)  =  (ONE  +  C2)*dnhat(i, 1) 
dnhat(i,2)  =  C2*dnhat (i , 1) 

nhat(i,3)  =  dnhat(i,2)  +  nhat(i,2) 

j_max  =  3 
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do  j  =4,  nfaces_j 


C2  =  max(  ONE,  (ONE  +  Cl*sin((j  -  2)  *  dpi))  ) 

!write(*,*)  i,  j,  C2 

dnhat(i, j-1)  =  min(  (C2*(nhat(i, j-1)  -  nhat(i, j-2))) ,  & 

((ONE  -  nhat(i, j-l))/(nfaces_j  -  j  +  D)  ) 

!write(*,*)  i,  j-1,  dnhat(i,j-l) 

nhat(i,j)  =  dnhat(i, j-1)  +  nhat(i,j-l) 

! write (*,*)  i,  j,  nhat(i,j) 

if  (  C2  >  C2_max  )  then 

C2_max  =  C2 
j_max  =  j 

end  if 

end  do 

if  (  nhat(i,nfaces_j)  /=  ONE  )  then 
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nhat(i,2:nfaces_j)  =  nhat(i,2:nfaces_j)  /  nhat (i ,nf aces_j ) 

end  if 

if  (  epO  /=  0  )  then 

do  j  =2,  nfaces_j 

ep  =  nhat(i,j)**2  *  (  1  -  nhat(i,j)  )  *  epO 

nhat(i,j)  =  (  ONE  -  ep  )  *  nhat(i,j)  +  fsh  *  ep 

! write (*,*)  'i’ , ’ j ^ , ’ep' , ’nhat' 

! write (*,*)  i,  j,  ep,  nhat(i,j) 

end  do 

end  if 

nxyf actor  =  noldf (i ,nf aces_j )  *  nshk(i) 

!nxyf actor  =  nshk(i) 

! nxyf actor  =  min(  nxyf actor,  C2_max  ) 

$  writedOl,  *)  'nxyf actor (i=l)  =  ’,  nxyfactor 
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do  j  =2,  nfaces_j 


nnewf(i,j)  =  nhat(i,j)  *  nxyf actor 
nnew(i,j-l)  =  HALF*(nnewf (i, j)  +  nnewf (i , j-1) ) 

end  do 

!  for  inflow  boundary  meeting  directly  with  body,  keeps  grid  from 
!  collapsing 

if  ( (algn_exclude  /=  0)  .and.  (i  <  algn_exclude) )  then 

nnewf (i,:)  =  nnewf (i+1,:) 
nnew(i,:)  =  nnew(i+l,:) 

end  if 

!  interpolate  solution  into  new  cell  center  positions 

call  interpextr (nold(i ,  1  :ncells_j )  jgrid’/oQ (R_ , i ,  1  :ncells_j )  ,nnew(i ,  1  :ncells_j )  ,xf] 
1 ,ncells_j , 1 ,ncells_j ) 

gridyoQ(R_ ,  i ,  1  :ncells_j )  =  xfnew(i ,  1  :ncells_j ) 

call  interpextr (nold(i ,  1  :ncells_j )  ,gridyoQ(RU_, i ,  1 : ncells_j )  ,nnew(i ,  1 : ncells_j )  ,x: 
l,ncells_j ,l,ncells_j) 

gridyoQ(RU_,  i ,  1  :ncells_j )  =  xfnew(i,  1  :ncells_j) 

call  interpextr (nold(i ,  1  :ncells_j )  jgrid’/oQ (RV_, i ,  1 : ncells_j )  ,nnew(i ,  1 : ncells_j )  ,x: 
1 ,ncells_j , 1 ,ncells_j ) 

gridyoQ(RV_,  i ,  1  :ncells_j )  =  xfnew(i,  1  :ncells_j) 
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call  interpextr (nold(i ,  1  :ncells_j )  jgrid'/oQ (E_ ,1,1  :ncells_j )  ,nnew(i ,  1  :ncells_j )  ,xf] 
1 ,ncells_j , 1 ,ncells_j) 

gridyoQ(E_ ,  i ,  1  :ncells_j )  =  xfnewd ,  1  :ncells_j ) 

call  interpextr (nold(i ,  1  :ncells_j )  ,gridyoQ(P_ ,1,1  :ncells_j )  ,nnew(i ,  1  :ncells_j )  ,xf] 
1 ,ncells_j , 1 ,ncells_j ) 

gridyoQ(P_ ,  1 , 1  :ncells_j )  =  xfnewd ,  1  :ncells_j ) 

call  interpextr  (noldd ,  1  :ncells_j )  ,gridyomud ,  1  :ncells_j )  ,nnewd ,  1  :ncells_j )  ,xfnei 
1 ,ncells_j , 1 ,ncells_j) 
gridymud ,  1  :ncells_j)  =  xfnewd,  1  :ncells_j) 

call  interpextr  (noldd ,  1  :ncells_j )  ,gridyolambdad ,  1  :ncells_j )  ,nnewd ,  1  :ncells_j ) ,: 
1 ,ncells_j , 1 ,ncells_j ) 

gridyiambdad ,  1  :ncells_j )  =  xfnewd ,  1  :ncells_j) 

gridyoQ(R_,CELL_MIN_INDEX:  (ncells_j+l)  ,CELL_M1N_INDEX:  (ncells_j+l)  )  & 

=  q_freestream(R_) 

gridy„Q(RU_,CELL_MIN_lNDEX: (ncells_j+l) ,CELL_MIN_1NDEX: (ncells. j+1) )  & 

=  q_freestreain(RU_) 

gridy„Q(RV_,CELL_MlN_lNDEX:  (ncells. j+1)  ,CELL_MIN_1NDEX:  (ncells. j  +  1)  )  & 

=  q.freestreain(RV.) 

gridyoQ(E.,CELL.MlN.INDEX:  (ncells. j+1)  ,CELL.M1N.1NDEX:  (ncells.j+1)  )  & 

=  q.freestreain(E.) 

gridy„Q(P.,CELL.MlN.INDEX:  (ncells.j  +  1)  ,CELL.M1N. INDEX:  (ncells.j+1))  & 

=  q.freestreain(P.) 

!  interpolate  face  center  positions... 

call  interpextr (noldf (i ,:) ,xf.c (i ,:) ,nnewf (i ,:) ,xfnew(i ,:) ,  & 

1 ,nf aces.j , 1 ,nf aces.j ) 

call  interpextr (noldf (i ,:) ,yf.c (i ,:) ,nnewf (i ,:) ,yfnew(i ,:) ,  & 

1 ,nf aces.j , 1 ,nf aces.j ) 


end  do 


if (vec_preserve  .eq.  1)  then 

do  i  =  1,  ncells_i 

!  Keep  all  points  on  same  vector  as  original  grid 

dxl  =  -gridy„A_i  (Y_,  i ,  1) 

dyl  =  grid°/oA_i(X_,i,l) 

dx2  =  xfnew(i,2)  -  xfnew(i,l) 

dy2  =  yfnew(i,2)  -  yfnew(i,l) 

na  =  (dxl*dx2  +  dyl*dy2)/((dxl*dxl  +  dyl*dyl) *sqrt (dx2*dx2 

! write (*,*)' dxl’ ,  dxl,  ’dx2’,  dx2,  ’dyl’,  dyl,  ’dy2’,  dy2, 

gridyox(i,2:nfaces_j)  =  na*nnewf (i,2:nfaces_j)*dxl  & 

+  grid'/oxCi,!) 

gridyoy(i,2:nfaces_j)  =  na*nnewf (i,2:nfaces_j)*dyl  & 

+  grid'/oyCi,!) 


end  do 

!  Keep  jmax  points  on  same  vector  as  original  grid 

dxl  =  -gridyoA_i(Y_,nfaces_i,  1) 

dyl  =  gridyA.i (X_ ,nf aces_i , 1) 

dx2  =  xfnew(ncells_i ,2)  -  xfnew(ncells_i , 1) 


+  dy2*dy2)) 

’na’ ,  na 
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dy2  =  yfnew(ncells_i ,2)  -  yfnew(ncells_i , 1) 

na  =  (dxl*dx2  +  dyl*dy2)/( (dxl*dxl  +  dyl*dyl) *sqrt (dx2*dx2  +  dy2*dy2)) 

! write (*,*) Mxl’ ,  dxl,  ’dx2',  dx2,  'dyl’,  dyl,  'dy2’,  dy2,  'na’ ,  na 

gridyox(nfaces_i,2:nfaces_j)  =  na*nnewf (ncells_i,2:nfaces_j)*dxl  & 

+  grid’/oxCnf  aces_i ,  1) 

gridyoy(nfaces_i,2:nfaces_j)  =  na*nnewf (ncells_i,2:nfaces_j)*dyl  & 

+  grid'/oyCnf  aces_i ,  1) 


else 

!  interpolate  interior  grid  points,  jmax  . . .  from  new  face  centers 

gridyox(2:  (nfaces_i-l)  ,2:nfaces_j)  & 

=  HALF* (xfnewd :ncells_i ,2 :nf aces_j )  +  xfnew(2:nfaces_i,2:nfaces_j)) 
gridyoy(2:  (nfaces_i-l)  ,2:nfaces_j)  & 

=  HALF* (yfnewd :ncells_i ,2 :nf aces_j )  +  yfnew(2:nfaces_i,2:nfaces_j)) 

!  extrapolate  x  coordinates  of  imin  surface,  y  coordinates  of  imax  surface 
!  (hardcoding  for  axisymmetric ,  positive  x  direction  flow,  orthogonal 
!  outflow  face . . . ) 

!  gridyoxd,2:nfaces_j)  =  TW0*xfnewd,2:nfaces_j)  -  gridyox(2,2:nfaces_j) 

!  gridyoy(nfaces_i,2:nfaces_j)  =  TW0*yfnew(ncells_i,2:nfaces_j)  & 

!  -  gridyoy(ncells_i,2:nfaces_j) 

!  Keep  imax  points  same  as  adjacent  to  keep  normal  to  flow. . . 
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gridy„x(l,2:nfaces_j)  =  gridyox(2,2:nfaces_j) 


!  Keep  jmax  points  on  same  vector  as  original  grid 

dxl  =  -gridyoA_i(Y_,nfaces_i,  1) 

dyl  =  grid’/oA.i  (X_  ,nf  aces_i ,  1) 

dx2  =  xfnew(ncells_i ,2)  -  xfnew(ncells_i , 1) 

dy2  =  yfnew(ncells_i ,2)  -  yfnew(ncells_i , 1) 

na  =  (dxl*dx2  +  dyl*dy2)/( (dxl*dxl  +  dyl*dyl) *sqrt (dx2*dx2  +  dy2*dy2)) 

! write (*,*)' dxl’ ,  dxl,  ’dx2’,  dx2,  ’dyl’,  dyl,  ’dy2’,  dy2,  ’na’ ,  na 

gridyox(nfaces_i,2:nfaces_j)  =  na*nnewf (ncells_i,2:nfaces_j)*dxl  & 

+  gridyox(nf  aces_i ,  1) 

gridyoy(nfaces_i,2:nfaces_j)  =  na*nnewf (ncells_i,2:nfaces_j)*dyl  & 

+  gridyoyCnf  aces_i ,  1) 


end  if 

Compute  the  grid  metrics 

do  j  =  1,  ncells_j 
do  i  =  1,  nfaces_i 

gridy„A_i(X_,i, j)  =  gridZyCi , j+1)  -  gridyoy(i,j) 
gridy„A_i(Y_,i, j)  =  -(gridy„x(i , j+1)  -  gridy„x(i, j)) 
gridy.A.i  (AMAG_ ,  i ,  j )  =  sqrt  (sum  (grid'/.A.i  ( 1 :  DOMAIN.DIM ,  i ,  j )  **2) ) 
end  do 
end  do 
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do  i  =  1,  ncells_i 
do  j  =  1,  nfaces_j 

gridy„A_ j  (X_ , i ,  j )  =  gridy„y(i,j)  -  gridy„y(i+l ,  j) 
gridy„A_ j  (Y_ , i ,  j )  =  -(gridy„x(i,  j)  -  grid'/oxCi+l ,  j) ) 
gridy.A.j  (AMAG_,i,  j)  =  sqrt (suni(gridy.A_j  (1  :DOMAIN_DIM, i ,  j) **2) ) 
end  do 
end  do 

do  j  =  1,  ncells_j 
do  i  =  1,  ncells_i 

dxl  =  gridy„x(i+l, j+1)  -  grid'/oxCi.j) 
dyl  =  grid'/oyCi+l, j+1)  -  grid'/oyCi,]) 
dx2  =  grid'/oxCi, j+1)  -  gridy„x(i+l, j) 
dy2  =  grid'/oyCi, j+1)  -  gridy„y(i+l, j) 
gridy.voKi, j)  =  HALF* (dxl *dy2  -  dx2*dyl) 
end  do 
end  do 

if (  Reynolds_num  /=  ZERO  )  then 
call  calc_viscous_metrics 0 
end  if 

!!$  writeClOl,  *)  'nshk(l)  =  nshk(l) 

!  !$ 

!!$  do  j  =  1,  ncells_j 

!!$  writedOl,  *)  'nnewf  (1 ,  d  j"''! d )  =  ’>  nnewf(l,  j+D 

!!$  writeClOl,  *)  'nnew(l,’>j>  ’)  =  d  nnew(l,j))  d^ho  =  grid'/oOCR.,  !>  j) 

!  !  $  end  do 
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! write (*,*)  ’Old  curvilinear  distance  from  wall  of  face  centers’ 
! write (*,*)  noldf ( : , : ) 

! write (*,*)  ’nhat  —  nondimensional  new  curvilinear  distance’ 

! write (*,*)  nhat(:,:) 

!write(*,*)  ’New  curvilinear  distance  from  wall  of  face  centers’ 
! write (*,*)  nnewf ( : , : ) 


!  !$ 
!  !$ 
!  !$ 
!  !$ 
!  !$ 
!  !$ 
!  !$ 
!  !$ 

!  !$ 


write(*,*) 
write(*,*) 
write (*,*) 
write (*,*) 
write (*,*) 
write (*,*) 
write (*,*) 
write(*,*) 


’Old  Face  Center  X-Coordinate ’ 
xf_c( : , : ) 

’Old  Face  Center  Y-Coordinate ’ 
yf_c(: , :) 

’New  Face  Center  X-Coordinate ’ 
xfnew( : , : ) 

’New  Face  Center  Y-Coordinate ’ 
yfnew( : , : ) 


close(lOl) 


deallocate (  xf_c  ) 
deallocate (  yf_c  ) 
deallocate (  xfnew  ) 
deallocate (  yfnew  ) 

!  deallocate (  x_c  ) 
!  deallocate (  y_c  ) 
deallocate (  dnx  ) 
deallocate (  dny  ) 
deallocate (  dnf  ) 
deallocate (  nold  ) 
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deallocate (  noldf  ) 
deallocate (  dnhat  ) 
deallocate (  nhat  ) 
deallocate (  nnew  ) 
deallocate (  nnewf  ) 
deallocate (  nlastcell  ) 
deallocate (  nlastface  ) 
deallocate (  nshk  ) 

end  subroutine  align_shock 

A.l  Shock  Location  Subroutine 

!  Aligns  Grid  to  and  clusters  points  around  shock,  keeping  adequate  boundary 

!  layer  resolution.  Taken  from  Gnoffo,  Hartung  and  Greendyke(AlAA-93-270-864) 


subroutine  shock_loc(s,  sshk,  comp) 

use  runtime_module 
use  constant_module 
use  ErrorOut_module 
implicit  none 


integer 

:  :  i,  j 

integer 

: :  ncells_i,  nfaces_i 

integer 

: :  ncells_j ,  nfaces_j 

integer 

:  :  cell_imin,  cell_imax 

integer 

:  :  cell_jmin,  cell_jmax 

!  counters 
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integer,  intent (in)  : :  comp  !  component  of  solution  vector  to  sense 

!  shock 

real (kind=RKIND) ,  dimension( :  ,  : ) ,  intent (in)  ::  s  !  cell  curvil.  dist.  array 
real (kind=RKIND) ,  dimension( : ) ,  intent(out)  ::  sshk  !  shock  location  array 

real (kind=RKIND)  ::  atest,  dtest,  propty,  proprty,  dtestref,  proprtyl 

ncells_i  =  num_i_cells (grid)  ! 

nfaces_i  =  num_i_f aces (grid)  ! 

cell.imin  =  CELL_M1N_INDEX  ! 

cell_imax  =  cell_max_i_index(grid)  ! 

ncells_j  =  num_j_cells(grid) 
nfaces_j  =  num_j_f aces (grid) 
cell.jmin  =  CELL_M1N_INDEX 
cell_jmax  =  cell_max_j_index(grid) 

sshk  =  ZERO 

select  case  (comp) 
case  (  R_  ) 

propty  =  ONE 
case  (  P_  ) 

propty  =  ONE/GAMMA 
case  default 

call  error2scr(" improper  shock  location  sense  variable:  must  be  p  or  rho") 


internal  cells  only 

faces  =  #  cells  +  n_ghost 

includes  ghost  cells  (0  or  1  -  nghost) 

includes  ghost  cells  (ncells_i  +  nghost) 
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end  select 


dtestref  =  fctrjmp  *  propty 

do  i  =  1,  ncells_i 

proprty  =  propty 

do  j  =  ncells_j ,  2,  -1 

proprtyl  =  proprty 
proprty  =  gridy„Q(comp,i,  j) 
dtest  =  dtestref  -  proprty 

if  (  dtest  <  -1 . 0e-6_RKIND  )  then 

atest  =  dtest  *  (  s(i,j+l)  -  s(i,j)  )  & 

/  (  proprtyl  -  proprty  )  +  s(i,j)  ! 

sshk(i)  =  atest  /  (  fsh  *  s(i,ncells_j)  )  ! 

! 

!write(*,  ’(A,  1,  A,  I,  A,  F8.3,  A,  F8.3)0  ' 

!write(*,  ’(A,  F8.3,  A,  F8.3,  A,  F8.3,  A,  F8. 

exit 

end  if 


interp  shk  s-value 

stretch/ compression 
factor 

s(’  ,  i,  ’ , ’ ,  j ,  0  =  ' ,  s(i,  j 
3)0  ’proprtyl  =  proprtyl, 
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end  do 


end  do 

! write (*,  ’(A)’)  'sshk  test  completed' 

! 

end  subroutine  shock_loc 

A. 2  Interpolation  and  Extrapolation  Subroutine 

!  Interpolation/extrapolation  subroutine  to  support  align_shock 
!  Arguments : 

!  xold,  yold,  xnew,  ynew  ...  x  indep  variable  y  dependent  var,  old  to  new 

!  [old  I  new] low  -  low  index  to  start  in  x. . .  [old | new] high 

!  assumes  x  is  increasing,  indices  are  increasing  from  low  to  high 

subroutine  interpextr (xold,  yold,  xnew,  ynew,  oldlow,  oldhigh,  newlow,  newhigh) 

implicit  none 

integer  ::  i,  j  !  counters 

integer,  intent (in)  : :  oldlow,  oldhigh,  newlow,  newhigh  !  indicial  extents 
integer  : :  oldlow2  !  temporary  lower  extent  of  old  data. . . . 

real (kind=RKlND)  : :  dxi  !  denominator  for  xold  difference 
real (kind=RKlND) ,  dimension( : ) ,  intent (in)  ::  xold,  yold,  xnew 
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real (kind=RKIND) ,  dimens i on( :) ,  intent (inout)  ::  ynew 


oldlow2  =  oldlow 

do  i  =  newlow,  newhigh 

if  (  xnew(i)  <=  xold(oldlow)  )  then  !  extrapolate  low  side 

dxi  =  ONE  /  (  xold(oldlow+l)  -  xold(oldlow)  ) 
ynew(i)  =  yold(oldlow)  & 

-  (yold(oldlow+l)  -  yold(oldlow))*dxi*(xold(oldlow)  -  xnew(i)) 

elseif  (xnew(i)  >=  xold(oldhigh)  )  then  !  extrapolate  high  side 

dxi  =  ONE  /  (  xold(oldhigh)  -  xold(oldhigh-l)  ) 
ynew(i)  =  yold(oldhigh)  & 

+  (yold(oldhigh)-yold(oldhigh-l) ) *dxi* (xnew(i)-xold(oldhigh) ) 

else  !  interpolate  from  low  to  high 

do  j  =  oldlow2,  oldhigh 

if  (  xnew(i)  <  xold(j+l)  )  then 

dxi  =  ONE  /  (  xold(j+l)  -  xold(j)  ) 
ynew(i)  =  yold(j)  & 

+  (yold(j+l)  -  yold(j ) ) *dxi* (xnew(i)  -  xold(j)) 


exit 
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else 


oldlow2  =  oldlow2  +  1 


end  if 


end  do 


end  if 


end  do 


end  subroutine  interpextr 
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